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L  INTRODUCTION. 

This  report  presents  some  navigation  and  almanac  algorithms  with  implementation  pro- 
grams that  are  written  in  the  Radio  Shack  TRS-so  MODEL  4  BASIC  language.  This  language 
can  be  customised  for  the  Navy  adopted  Sharp  PC-1600A  or  other  BASIC  language  com- 
puters. 

Because  the  PC-1600A  has  limited  memory,  many  of  the  classical  navigation  formulas 
from  spherical  trigonometry  have  been  rewritten  in  a  form  which  minimises  or  eliminates 
the  need  to  determine  special  cases— especially  those  of  quadrant.  Details  are  contained 
in  Section  H 

The  navigation  algorithms  are  detailed  in  the  BASIC  program  NAVALGOR,  presented  in 
Section  IE.  Programs  for  the  determination  of  the  position  of  the  Sun,  Moon,  Venus,  Mars, 
Jupiter  and  Saturn  are  contained  in  the  BASIC  program  NAVBPHM  which  is  presented  in 
Section  IV.  These  programs  reproduce  the  tables  in  The  Nautical  Almanac  and  The  Air 
Almanac  to  within  0.2*. 

IL  THE  GENERAL  SPHERICAL  TRIANGLE 
AND  COORDINATE  TRANSFORMATIONS. 

A*  Background.  Many  of  the  algorithms  and  procedures  used  for  navigational  compu- 
tations were  designed  for  use  with  tables  of  sines,  cosines  and  tangents.  To  save  space  and 
eliminate  repetition,  these  tables  usually  contained  only  values  for  the  first  quadrant.  To 
generate  values  for  other  quadrants,  many  cumbersome  rules  were  adopted.  These  rules — 
contained  in  classical  references  such  as  Bowditch  [Ref.  1] — are  promulgated  in  modern 
computer  programs,  causing  needlessly  additional  programming  effort.  This  additional  ef- 
fort can  largely  be  eliminated  by  using  the  concept  of  the  general  spherical  triangle  which 
is  described  in  the  next  section. 

B.  Theory.  The  formulas  for  the  general  spherical  triangle  expounded  on  by  Wm.  Chau- 
venet  [Ref.  3)  in  his  book,  A  Treatise  on  Plane  and  Spherical  Trigonometry,  which  was  first 
published  in  1850.  One  usually  thinks  of  a  spherical  triangle  as  looking  something  like  the 
object  depicted  in  Figure  1,  where  the  arc  length  of  each  leg  is  less  than  180°.  Chauvenet  's 
general  triangle  does  not  have  the  180°  limitation  and  so  objects  similar  to  the  one  shown 
in  Figure  2  are  also  considered  to  be  spherical  triangles.  In  Figure  2,  the  crosspoint  P 
is  180°  from  the  vertex  of  angle  A.  In  labeling  spherical  triangles  it  is  important  that  an 


Figure  1.  A  Spherical  Triangle 


Figure  2.  A  General  Spherical  Triangle 

angle  be  labeled  with  a  capital  letter  and  that  the  opposite  side  be  labeled  with  the  same, 
but  lower  case  letter. 

Depending  upon  the  quantities  to  be  determined,  either  the  set  of  equations 


cos  a  =  cos  b  cos  c  +  sin  b  sin  c  cos  A, 
sin  a  cos  B  -  cos  6  sin  c  -  sin  b  cos  c  cos  A, 
sin  a  sin  B  =  sin  b  sin  A, 


(i) 


or  the  polar  form  of  Equs.  (1) 


cosA  =  -  cos  B  cos  C  +  sin  B  sin  C  cos  a, 
sin  A  cos  6  =  cos  5  sin  C  +  sin  B  cos  C  cos  a, 
sin  A  sin  6  =  sin  B  sin  a, 
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(2) 


may  be  used. 

Chauvenet  stated  that,  "It  will  be  found  that  all  the  six  cases  of  the  general  triangle 
admit  two  solutions,  but  that  then  all  become  determinate,  when,  in  addition  to  the  other 
data,  the  sign  of  the  sine  or  cosine  of  one  of  the  parts  is  given.* 

Formulas  for  spherical  coordinates  should  be  developed  so  that  quantities  such  as 
latitude,  declination  and  altitude,  which  have  values  in  the  range  -90°  to  +90°,  are  de- 
termined by  arcsine  or  arctangent  formulas;  so  that  quantities  such  as  longitude,  asimuth, 
Greenwich  hour  angle  and  local  hour  angle,  which  have  values  in  the  range  -180°  to  +180° 
or  0°  to  360°,  are  determined  by  both  a  sine  and  a  cosine  formula  used  in  conjunction  with 
a  quadrant  determining  arctangent  function  (the  qatn  function);  and  so  that  a  quantity 
such  as  distance,  which  has  a  spherical  arc  value  in  the  range  of  0°  to  180°  (assuming  the 
user  wishes  the  shortest  distance),  is  determined  by  an  arccosine  formula. 

Illustration  1.  Consider  the  "inverse9  problem  of  spherical  geometry,  which  is:  Given 
the  latitude,  ^,  and  longitude,  A,  of  two  points  Pi{fa,Xi)  and  JV^a,  A3),  determine  the 
distance  d,  the  forward  asimuth  ct\2  and  the  backward  asimuth  a?i.  The  solutions  will  be 
derived  using  the  convention  that  eastern  longitudes  are  positive  and  western  longitudes 
are  negative.  Also,  northern  latitudes  are  positive  and  southern  latitudes  are  negative. 
The  geometry  is  illustrated  in  Figure  3. 

North  Pole 


PiOiAi)      ^^-—_  ±^    PiihM 


Figure  3.  Geometry  of  the  Navigation  Triangle. 
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Solve  for  d  and  Q13.  Label  the  angles  and  legs  so  that  A  =  A?  -  Ai,  B  =  crn, 
(7  =  360°  -  a31,  a  =  d,  b  =  90*  -  fa  and  c  =  90°  -  fa .  Substituting  these  into  Equs.  (1) 
and  simplifying  the  differences  of  angles,  we  obtain 

cos d  =  sin  fa  sin fa  +  cos^2Cos^icos(A3  -  Ai), 
sinicosais  =  sin  fa  cos  fa  —  cos^sin^iCosfA?  —  Ai), 
sinisinaia  =  cos^2sin(As  -  Ai). 

The  first  equation  can  be  used  to  solve  for  d.  Using  the  principle  angle  of  the  arccosine,  d 
lies  in  the  range  0°  to  180°.  If  the  non-principle  angle  is  used,  d  lies  in  the  range  180°  to 
360°,  we  can  either  travel  the  shortest  great  circle  distance  from  Pi  to  Pi  or  we  can  travel 
the  great  circle  route  which  goes  around  the  backside  of  the  earth — the  former  solution  is 
assumed  to  be  preferred.  Hence,  d  is  restricted  to  lie  between  0°  and  180°.  Multiply  d  by 
60  to  get  the  distance  in  nautical  miles. 

Chauvenet  states  that  this  system  of  equations  has  two  solutions,  but  has  only  one 
solution  if  the  sign  of  the  sine  or  cosine  of  one  of  the  parts  is  known.  Since  d  lies  between 
0°  and  180°,  the  sign  of  sin  d  is  known,  thus  the  system  has  only  one  solution,  further, 
and  most  important,  the  sign  of  sin  d  m  positive  for  all  d  between  0°  and  180°  and  so  there 
are  no  special  cases. 

The  first  equation  is  used  to  compute  d: 

d  =  9iccos[8mfa8mfa+cosfaco6faco6(X^-Xi)\.  (3) 

The  second  and  third  equations  are  used  to  compute  an.  A  natural  tendency  would  be  to 
solve  either  the  second  equation  or  the  third  equation,  but  doing  so  loses  the  information 
concerning  the  quadrant  of  an .  Use  the  second  equation  to  determine  cos  an  and  the  third 
equation  to  determine  sin  an  ■  With  both  the  sine  and  the  cosine  known,  the  quadrant  is 
uniquely  determined  using  the  quadrant  determining  arctangent  function.  The  solution  is 

an  =  qatn(sin  an,  cos  cm) 

as  long  as  d  is  not  0°  or  360°.  Even  though  the  chance  of  division  by  zero  with  the 
occurrence  of  a  value  for  d  of  exactly  0°  or  exactly  360°  is  very  remote,  it  can  be  eliminated 
entirely.  Since  the  sign  of  sin  d  is  positive  or  zero  in  both  the  second  and  third  equations, 
it  effectively  cancels  out  using  the  qatn  function.  The  asimuth  a\2  is  found  solving 

<*ia  =qatn[cosfosin(A3-Ai),sin<fecos^i  -cos ^3 sin ^i  cos(A3  -  Ai)|.  (4) 
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To  find  the  back  azimuth,  label  the  angles  and  legs  so  that  A  =  A?-  Ai ,  B  =  360°  -arsi , 
C  =  cti2,  a  =  <i,  6  =  90°  ~<£i  and  e  =  90°  -  fo.  Substitute  these  into  Equs.  (1)  and  simplify 
to  obtain 


<*3i  =qatn[-cos^isin(Ai  -  As),sin^icos<fa  -cos<h  sin<focos(Ai  -As)]. 


(5) 


Hhistration  2*  Consider  the  'direct*  problem  of  spherical  geometry,  which  is:  Given 
the  latitude,  <j>,  and  longitude,  A,  of  a  point  <Pi(^i,Ai),  the  distance  d  and  the  forward 
azimuth  ai3  to  a  point  /s(<fe,  As),  determine  <fo,  As  and  the  backward  azimuth  a?i. 

Solve  for  <fo  and  As.  Relabel  the  angles  and  legB  so  that  A  =  oris,  5  =  Ai  -  As, 
C  =  360°  -  atis,  a  =  90°  -  fa,  b  =  d  and  c  =  90°  -  ^.  Substituting  these  into  Equs.  (1) 
and  simplifying  the  differences  of  angles,  we  obtain 

sin ^s  =  cosdsin^i  -f  sindcos^i  coscria, 
cos^scos(A3  -  Ai)  =  coadcosfa  -sindsin^icosais, 
cos  fa  sin(As  -  Ai )  =  sin  isin  oris. 

The  first  equation  can  be  used  to  solve  for  fa.  Using  the  principle  angle  of  the  arcsine,  fa 
lies  in  the  range  —90°  to  90°,  which  is  correct  for  a  latitude.  Since  fa  lies  between  -90° 
and  +90°,  the  sign  of  cos  fa  is  known,  thus  the  system  has  only  one  solution.  Further,  and 
again  most  important,  the  sign  of  cos  fa  is  positive  for  all  fa  between  —90°  and  +90°  and 
so  there  are  no  special  cases. 

The  first  equation  is  used  to  compute  fa: 

fa  =  arcsin(cos<iain  fa  +  awdco&fa  coscris).  (6) 

The  second  and  third  equations  determine  cos(As  -  Ai)  and  sin(A3  -  Ai).  As  before,  with 
both  the  sine  and  the  cosine  known,  the  quadrant  is  uniquely  determined  using  the  qatn 
function.  Since  cos  fa  is  always  positive,  it  can  be  eliminated  in  both  equations.  Then 

As  =  Ai  +qatn(sin<i8in  oris,  cos  <f  cos  ^i  -sinrfsin^icosais).  (7) 

The  back  azimuth  can  be  found  by  relabeling  the  angles  and  legs  so  that  A  =  asi  , 
B  =  360°  -  oai,  C  =  Ai  -  As,  a  =  90°  -  ^s,  6  =  90°  -  ^i  and  c  =  d.  Substituting  these 
into  Equs.  (1)  and  simplifying  the  differences  of  angles,  we  obtain 

&2i  =  qatn(-  cos  fa  sin  a i3,sin<^i  sin  d-  cos  fa  cos d  cos  ais).  (8) 
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This  second  formula  for  the  back  azimuth  is  used  in  the  spirit  that  the  unknown  parameters 
should  be  determined  as  a  function  of  the  known  parameters  rather  than  determining 
the  second  unknown  parameter  as  a  function  of  the  previously  evaluated  first  unknown 
parameter.  That  is,  any  one  of  the  unknowns  can  be  found  as  a  function  of  the  knowns 
without  having  to  first  find  one  of  the  other  unknowns. 

Hhutration  3.  Although  derived  from  different  principles,  the  equations  of  coor- 
dinate transformation  obey  the  rule  of  Chauvenet.  Fbr  example,  the  equations  which 
transform  local  hour  angle  and  declination  to  altitude  and  aiimuth  are 

cosasinA  =  —  cos  £  sin  A, 

cos  a  cos  A  =  sin  6  cos  ^  -  cos  6  sin  ^  cos  h,  (9) 

sin  a  =  sin  6  sin  $  +  cos  6  cos  $  cos  A, 

where  a  =  altitude,  A  =  aiimuth,  6  =  declination,  h=  local  hour  angle  and  <j>  =  latitude 
[Ret  4,  pg.  26]. 

The  third  equation  is  used  to  solve  for  altitude 

a  =  arcsin(sin$sin^  +  cos$cos^cosfc).  (10) 

The  altitude  lies  between  -90°  and  +90°,  and  so  is  correctly  determined  by  the  arcsin 
function.  Since  cos  a  is  always  positive  it  can  be  eliminated  between  the  first  two  equations 
when  using  the  qatn  function.  The  local  hour  angle,  which  lies  in  the  range  0°  to  360°,  is 
determined  by 

A  =  qatn(-cc«$smA,sin$cos^-cos$sin^cos/i).  (11) 

C.  Simplification  of  Standard  Works.  Using  the  general  triangle  and  the  qatn  func- 
tion, the  simplifications  listed  below  can  be  made  in  Volume  II  of  Bowditch  [Ref.  1].  The 
various  sections  are  preceded  by  the  section  symbol  §  and  equation  numbers  from  Bowditch 
within  those  sections  are  prefixed  with  the  letter  B.  Unprefixed  equation  numbers  are  those 
contained  in  this  document. 

§700.  Solving  for  altitude. — Using  the  conventions  that  northern  latitudes  and 
declinations  are  positive  and  that  southern  latitudes  and  declinations  are  negative,  the 
altitude  of  a  star  can  be  computed  directly  from  Equ.  B(2a)  or  Equ.  B(2b).  Equ.  B(2a)  is 
equivalent  to  Equ.  (10).  All  of  the  special  cases  are  thus  eliminated. 

§707.  Solving  for  aiimuth.— Instead  of  using  Equ.  B(4b)  or  Equ.  B(5b),  use 
Equ.  (11).  See  §708  below. 
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§703.  Tim©  aiimuth* — Bowditch  uses  meridian  angles  which  can  be  labeled  either 
east  or  west.  Since  altitude/ aiimuth  Is  a  left-handed  coordinate  system,  the  convention  is 
that  west  meridian  angles  are  positive  and  east  meridian  angles  are  negative. 

Example  L— The  latitude  of  the  observer  is  30°25/0N;  the  declination  of  the  celestial 
body  is  22°06.'2N;  the  local  hoar  angle  is  39°54.'7W.  Using  Equ.  (11), 

A  m  qatn(-co822°06.'2sin  39°  54.'7,sin  22d06.'2cos30°25.'0 

-  cos22o06.'2sin30o25.'0cos39°54.'7) 

=  qatn[-(0.92651)(0.64161),  (0.37628)  (0.86237) 

-  (0.92651)(0.50628)(0.76703)] 
=  qatn(-0.59445, 0.32449  -  0.35980) 
=  qatn(-0.59445,  -0.03531) 

=  -93.°39913 
=  266.°60087 
A  =  266°36.'l. 

Example  2— The  latitude  of  the  observer  is  30°25.'0S;  the  declination  of  the  celestial 
body  is  22°06.'2N;  the  local  hour  angle  is  39°54.'7E.  Using  Equ.  (11), 

A  =  qatn(-cos22o06.'2sin(-39o54.'7),sin22o06./2cos(-30o25/0) 

-  cos22°06.'2sin(-30o25/0)  cos(-39°54./7)] 
=  qatn[-(0.92651)(-0.64161),  (0.37628)  (0.86237) 

-  (0.92651)(-0.50628)(0.76703)] 
=  qatn[0.59445, 0.32449  -  (-0.35980)] 
=  qatn(0.59445, 0.68429) 

=  40.°98141 
A  =  40°  58/9. 

D.  Spheroid  Earth  and  Great  Circle  Formula*. 

1.  Spheroid  Earth  Formulas.  The  spheroid  earth  algorithms  are  those  of  Thomas 
[Ref.  5].  Some  changes  by  Shudde  [Ref.  6]  make  quadrant  determination  automatic  by 
use  of  the  qatn  function.  East  longitudes  and  south  latitudes  are  negative.  The  earth's 
equatorial  radius  is  denoted  by  a9  and  the  earth's  flattening  factor  is  denoted  by  /. 


a.  Inverse  Solution.  Input  variables  are  the  latitude  fa  and  longitude  Ai  of  Pi, 
the  latitude  fa  and  longitude  Aa  of  Pa,  fa  and  A3.  Output  variables  are  distance,  (S/a*), 
forward  asimuth  from  Pi  to  P3,  cria,  and  back  asimuth  from  P2  to  Pi,  otai-  All  angular 
input  and  output  is  in  radians.  Compute: 

$i  =  tan"1  [(1  -  /)  tanfc],  for  t  =  1, 2, 
AA  -  Aa  -  Ai,  9m  =  (fa  +  fa)/2,  A*m  =  (fa  -  fa) /*, 
S  m  cos3  A^m  -  sin3  §m  L  m  sin3  Atfm  +  JJ  sin3(AA/2), 
<*  =  2  sin"1  (I1/3),  U  =  2sin3Cco83  A0m/(1  - X), 
y  =  2sm3A^cos30m/I,  JT=CT  +  7,  7  =  (7-7, 
r  =  <J/sin(i,  0  =  4T3,  £  =  2cosci,  ;4  =  0£, 

c  =  r-(i-fl/2,  m=x{A+cx), 

B  =  20,  n,  =  F (B  +  £Y),  n3  =  DXY,  (12) 

M  =  f{TX  -  r)/4,  M  =  f{nx  -  na  +  n3)/64, 
5/a,  =  (T  -  M  +  fa<i) sine*,  M  =  32T  -  (20T  -  A) X  -  (B  +  4)7, 

P  =  2K  -  £(4  -  X),  G  =  /r/2  +  /^ilf/H  Q  m  -(PG tan  AA)/4, 
AA'm  =  (AA+Q)/2, 

1 1  =  qatn(-  sin  A0m  cos  AA Jn,  cos 0m  sin  AA^), 
*a  =  qatn(cos  A0m  cos  AA^ ,  sin  $m  sin  A  AJJ , 
cfia  =  ti  +  fa  and  aai  =  ii  — 1%. 

b.  Direct  Solution.  Input  variables  are  the  latitude  fa  and  the  longitude  Ai  of 
Pi,  the  forward  asimuth  aia  from  Pi  to  Pa,  and  the  distance  (S/ae)  from  Pi  to  Pj.  The 
output  variables  are  the  latitude  fa  and  longitude  A?  of  Pa,  and  the  backward  asimuth 
arai  from  P3  to  Pi.  All  angular  input  and  output  is  in  radians.  Compute: 

fa  =  tan"1  [(1-/)  tan*], 
M  =  cosfa  sinotia,  c\  =  /Af,  Ca  =  /(l  -  M2)/*, 
D  =  (l-  ca)(l  -  ca  -  cjAf),  P  =  ca(l  +  ciM/2)/Dt 
N  =  cosfa  cos  aria,  <T\  =  qatn(#,sinfa), 
d  =  (S/a9)/D,  ti  =  2((Ti  -  d),  W  =  1  -  2Pcost», 
7  =  cos(u  +  <J),  X  =  c£  sin  icos(i(2K3  -  1), 
r  =  2PVWsin<i,  A*  =  d  +  X  -  K,  (13) 
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K  m  [M*  +  (tfcos  Aa  -  sin*!  sin  A*)a]1/3, 
tan  $2  =  (sin  0i  cos  A<r  +  iV  sin  Atr)/Kj 
^  =  tan-1[(tan^)/(l-/)], 

Aj?  as  qatn(sin  Ac  sin  aia ,  cos  $i  cos  A<r  -  sin  0i  sin  A<r  cos  aia), 
IT  =  Ci(l  -  ca)A<T  -  CiCa sin  A<tcos(2(7i  -  A<r), 
X7  m  Xi  +  Aij  -  iT, 
a^i  =  qatn(-M,  -(N  cos  A(7  -  sin  9\  sin  A<r)]. 

2.  Spherical  Earth  Formulas,  The  inverse  solution  [Equs.  (3),  (4)  and  (5)]  and  direct 
solution  [Equs.  (6),  (7)  and  (8)]  formulas  were  developed  in  Section  B  of  this  report.  They 
are  summarised  here  for  convenience. 

a.  Inverse  Solution.  Input  variables  are  the  latitude  fa  and  longitude  \\  of  Pi ,  the 
latitude  fa  and  longitude  A3  of  P3,  fa  and  A3.  Output  variables  are  distance,  d,  forward 
asimuth  from  Pi  to  Pj,  a  13,  and  bach  asimuth  from  Pj  to  Pi,  0:31.  All  angular  input  and 
output  is  in  radians.  Compute: 

d  =  arccosfsin  fa  sin  fa  +  cos  fa  cos  ^1  cos(A?  -  Ai )], 
aia  =  qatn(cos^jsin(Aa  -  Ai),sin<^jcos^i -cos^sin^icos(Aa  -  Ai)),  and  (14) 

otai  =qatn(-cos^isin(Ai  -  A3),  sin  ^1  cos  ^3  -cos^isin^3Cos(Ai  -A3)). 

b.  Direct  Solution.  Input  variables  are  the  latitude  fa  and  the  longitude  Ai  of  Pi, 
the  forward  asimuth  a\i  from  Pi  to  P3,  and  the  distance  d  from  Pi  to  P3.  The  output 
variables  are  the  latitude  fa  and  longitude  A3  of  P3,  and  the  backward  asimuth  0:31  from 
P3  to  Pi .  All  angular  input  and  output  is  in  radians.  Compute: 

fa  =  arcsin(cosdsin^i  +  sin  d  cos  ^1  cos  0:13), 

A3  =  Ai  +qatn(sin<isin  0:13,  cos  d cos  <0i  -  sin  d  sin  ^1  cos  an),  and  (15) 

<*2i  =  qatn(—  cos  ^1  sin  a12,  sin  fa  sin  d  -  cos  fa  cos  d  cos  an). 

c.  Given  Longitude,  Find  Latitude.  Input  variables  are  the  latitude  fa  and 
longitude  Ai  of  Pi,  the  forward  asimuth  0:13  from  Pi  to  P3,  and  the  longitude  A  of  some 
point  P  on  the  great  circle  joining  Pi  and  P3.  The  output  variable  is  <j>,  the  latitude  of 
P.  All  angular  input  and  output  is  in  radians.  Using  Equs.  (2),  label  the  angles  and  legs 
so  that  A  =  360°  -  <*2lj  B  =  0:13,  C  =  A  -  Ai,  a  =  90°  -  fa,  b  =  90°  -  ^  and  c  =  d. 
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Substituting  these  into  Equs.  (2)  and  amplifying  the  differences  of  angles,  we  obtain 

coscKai  =  -cosai2Cos(A-Ai)  +  sinai2  8in(A- Ai)sin^i, 
sin crai sin^  =  cosai2sin(A  -  Ai)  +  sin atu cos(A  -  Ai)sin^i,  and 
sina^i  cos^  =  sin  aria  cos  <^i .. 

Since  the  only  quantity  of  interest  is  ^,  and  since  latitude  lies  in  the  range  of  -90°  to  +90° , 
it  suffices  to  divide  the  second  equation  by  the  third  equation  to  obtain  an  expression  for 
tan  <f>>  We  obtain 

"cos  oris  sin(A~Ai)+sincri2Cos(A- Ai)sin^i 


<£  =  tan 


-1 


sin  ai2  cos  <j>\ 


(16) 


For  almost  all  applications,  Equ.  (16)  will  suffice.  It  is,  however,  possible  for  the  denom- 
inator of  Equ.  (16)  to  become  zero  in  certain  circumstances.  We  now  face  a  situation  in 
which  special  cases  must  be  examined. 

One  way  to  handle  the  dilemma  is  to  first  evaluate  the  denominator  and  the  numerator. 
If  the  value  of  the  denominator  is  zero,  then  ^  =  90°  if  the  numerator  is  positive  or 
$  =  -90°  if  the  numerator  is  negative.  If  the  denominator  is  not  aero,  then  use  Equ.  (16). 

A  second  way  to  handle  the  dilemma  is  to  compute 

^  =  qatn[cos ayx sin(A  -  Ai)  +  sin oti2 cos(A  -  Ai)  sin ^1, sin <x\i  cos ^1],         (17) 

but  this  leads  to  another  problem  because  ^  may  now  be  in  the  range  of  - 180°  to  + 180° .  To 
correct  the  output  to  a  proper  latitude,  subtract  180°  if  0  >  90°  or  add  180°  if  ^  <  -90°. 

The  question  of  which  method  to  use  is  a  matter  of  personal  choice.  Generally,  a 
simple  one  line  logic  function  can  handle  all  of  the  special  cases. 

d.  Find  the  Vertex*  The  vertex  of  a  great  circle  route  is  the  most  northerly  or/and 
the  most  southerly  point  on  the  great  circle.  Bowditch  [Ref.  1,  §1016,  Equs.  B25-B30]  gives 
formulas  for  first,  computing  the  latitude  of  the  vertex,  and  second,  computing  the  longi- 
tude of  the  vertex  as  a  function  of  the  latitude  of  the  vertex.  The  Bowditch  latitude  formula, 
from  a  straight  forward  application  of  the  Law  of  Sines,  is  <^,  =  cos-1  (sin  0:12  cos  <t>\ ),  where 
<j>*  is  the  latitude  of  the  vertex  and  <pi  and  0*13  have  the  same  meaning  as  in  the  previous 
section.  The  principle  angle  of  the  arccosine  lies  between  0°  and  180°  whereas  latitudes 
lie  between  -90°  and  +90°.  Bowditch  illustrates  only  cases  in  which  <j>*  is  between  0°  and 
+90°  and,  unfortunately,  makes  no  mention  of  how  to  handle  cases  when  the  arccosine  is 
greater  than  90°. 
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There  are  several  alternate  methods  of  developing  equations  for  the  vertices.  One 
method  is  presented  here  and  an  alternate  method  is  presented  at  the  end  of  §•.  Input 
variables  are  the  latitude  $\  and  longitude  Ai  of  Pi,  and  the  forward  asimuth  a^  from  Pi 
to  Pa.  Let  the  latitude  and  longitude  of  some  arbitrary  point  P  on  the  great  circle  from 
Pi  to  Pa  be  denoted  by  <j>  and  A,  respectively.  All  angular  input  and  output  is  in  radians. 
Label  the  angles  snd  legB  so  that  A  =  A  -  Ai,  B  =  ati3,  C  =  360°  -  otqi  ,  a  =  d,b  =  90°-^ 
and  c  =  90°  -  ^i.  Substituting  these  into  Equs.  (1)  and  simplifying  the  differences  of 
angles,  we  obtain 

cos  d  =  sin  j>  sin  <£i  +  cos^cos^icos(A-  Ai), 
sindcosaia  =  sin^cos^i  ~cos^ain^icos(A  -  Ai),  and 
sin  <i  sin  an  =  cos^sin(A  -  Ai). 

Divide  the  third  equation  by  the  second  equation  to  obtain 

cos^sin(A-  Ai) 


tan  a\2  =  — 


sin^cos^i  —  cos^sin^ico8(A  —  Ai)* 
Then  rearrange  the  equation  in  the  form 

tanaufsin^cos^i  -cos^sin^icos(A-  Ai)]  =  cos^sin(A  -  Ai).  (18) 

Next,  we  wish  to  find  what  value  of  ^  is  an  extremum  of  A.   To  do  this  we  must  first 
compute  d<t>/d A  in  Equ.  (18): 


tan  aria 


cos^cos^i^ry+sin^sin^i  cos(A  -  Ai)^y  +  cos^sin^i  sin(A  -  Ai) 


=  -sin^8in(A  —  Ai)^y  +  cos^cos(A  —  Ai). 


Then  we  must  set 


Thus, 


dx 


tanai2Cos^vsin^isin(Av  -  Ai)  =  cos^,  cos(A«  -  Ai), 
which  rearranges  to 


tan(A,  -  Ai)  =  -: 


1 


or 


A,  =  Ai  +  tan"1 


sin  <j>\  tan  a\q 

1 
sin^i  tanau 


(19) 
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Once  A,  is  known,  fa  can  be  found  using  one  of  the  methods  in  §c  such  as  Equ.  (17).  Since 
|A»  —Ai|  <  90°,  it  follows  that  fa  must  be  in  the  same  hemisphere  as  fa.  Thus,  if  ^i  is  in 
the  northern  (southern)  hemisphere,  then  fa  must  be  the  northern  (southern)  vertex. 

e.  Given  Latitude,  Find  Longitude.  Input  variables  are  the  latitude  fa  and 
longitude  Ai  of  Pi,  the  forward  asimuth  an  from  Pi  to  P3  and  some  latitude  fa  It  is 
convenient  but  not  necessary  to  know  the  latitude  of  the  vertex  fa  of  the  great  circle  from 
Pi  to  Pj.  We  wish  to  find  the  longitude(s)  on  the  great  circle  from  Pi  to  Pj  which  have 
latitude  fa  There  are  several  possible  outcomes:  (1)  if  \<j>\  >  |^,|,  there  is  no  solution;  (2) 
if  <p  -  fa%  there  is  one  solution  and  that  is  A  =  A,;  and  (3)  if  |^|  <  |^,|,  there  are  two 
longitudes  which  satisfy  the  given  conditions. 

Rewrite  Equ.  (16)  in  the  form 

cos  aia  cos  <j>  sin(A  -  Ai )  +  sin  oti3  sin  fa  cos  ^  cos(  A  -  Ai )  =  sin  ai3  cos  fa  sin  fa 

or 

5sin(A  -  Ai)  +  Ccos(A  -  Aj)  =  K,  (20) 

where  S  =  cosaijcos^,  C  =  sin 0:12 sin ^1  cos^,  and  K  =  sin  am  cos^i  sin^.  Define  p  >  0 
and  rf  as  the  solutions  to  the  equations 


C  =  pcoarf    and    S  =  painty.  (21) 


Solving,  we  find  that 


/>  =  +V$2  +  C2    and 
ff  =  qatn($,C) 
Substituting  Equs.  (21)  into  Equ.  (20)  and  rearranging,  we  find 

cos(A  -  (Ai  +  ri)\  =  K/p. 

So  that 

A-(Ai+i?)  =  ±cos-1(A'//?), 

or 

A  =  (Ai  +  ri)  ±  cos"1  (K/p) .  (22) 

Recall  that  a  knowledge  of  fa  is  not  known  to  determine  the  number  of  solutions.  (1)  if 
\K/p\  >  1,  there  are  no  solutions;  (2)  if  \K/p\  =  1,  there  is  one  solution;  and  (3)  if 
\K/p\  <  1,  there  are  two  solutions. 
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In  the  previous  section  it  was  mentioned  that  there  was  yet  another  way  to  determine 
the  position  of  the  vertex.  To  do  this,  for  fixed  <£,  define  /(A)  in  terms  of  Equ.  (20)  as 

/(A)  =  5sin(A  -  Ai)  +  Ccos(A  -  Ai)  -  K. 

Then  /'(A),  the  derivative  of  /(A)  with  respect  to  A  is 

f(X)  =  5cos(A  -  At)  -  Csin(A  -  \t) 

Define  A*  so  that  /*(A*)  =  0.  Thus 

5cos(A*  -  At)  -  Csin(A*  -  At)  =  0, 

OT  sin(A*-Ai)  _  S 

cos(A*-Ai)      C" 
or 

A*-Ai  =  qatn(5,C), 

so  that 

A*  =  Ai  +  qatn(S,  C)  =  Ai  +  tj. 

Hence  A*  =  Ai  +  rj  is  an  extremnm  of  Equ.  (20),  or  A,  =  A*.  Note  that 

/"(A*)  =  -Ssin(A*  -  A0  -  Ccos(A*  -  Ai) 
a  —  Ssiatf  —  C  cos  if 

p      p 

S2  +  C2 
P 
=  -p  <  0. 

Thus  A,  =  A*  is  the  northern  vertex. 

E.  Rhumb  Line  (Mercator)  Formulas. 

1.  Spheroid  Earth  Formulas.  Bowditch  [Ref.  1,  Explanation  of  Table  5,  Meridional 
Parts]  gives  the  formula  for  the  computation  of  meridional  parts  as 

M s  a.  In  tan  f  45°  +  ^  -a,  (Vsin^  +  jsin3^+  jsin5^  +  ...] 

=  aelntan(45°  +  |)  -a*e  (esin<£  +  -r-sin3^  +  ■=- ain5<^+  •••) 
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where  M  is  the  number  of  meridional  parts  between  the  equator  and  the  given  latitude 
4,  a,  is  the  equatorial  radius  of  the  earth,  and  t  =  y/(2f  -  P)  is  the  eccentricity  of  the 
earth,  where  /  =  1/298.26  is  the  WGS  1972  earth  flattening  factor. 
Using  the  identity  [Ref.  7,  #769) 

*(i^)=l[*+i*J+k+?*T+"]'  i*J<1i 

it  is  possible  to  rewrite  the  formula  for  M  as 

M-*h*(«.+i)-¥*(£3$) 


=  a« 


«/3' 


=  aeln 


/l  +  esinAg/a" 
^l-esinjy 


If  we  denote  the  course  (forward  asimuth)  by  a13,  then  [Bowditch,  Vol  EL,  §1001  and 
§1013]  AA  =  Aq  -  Ai  =  mtanau  where  m  =  M(fo)  -  M(^i),  so  that 

(Aq  -  Ai)/Oe 


tan  <x\2  = 


fc -(**»)  .h^£B 

/l  +  esinfoV7*         /l  +  gsin^y^ 
yl-esin^y  yl-esin^iy 


(23) 


This  equation  is  also  given  by  Shufeld  and  Newcomer  [Ref.  8,  pp  81-84).  It  is  the  spheroid 
earth  rhumb  line  equation. 

a.  Inverse  Solution.  Input  variables  are  the  latitude  <fi\  and  longitude  Ai  of  Pi ,  the 
latitude  <j>i  and  longitude  Aq  and  P^}  fa  and  Aq.  Output  variables  are  distance,  d,  forward 
asimuth  aria  from  Pi  to  Pj,  and  back  asimuth  an  from.  Pj  to  Pi .  All  angular  input  and 
output  is  in  radians.  In  Equ.  (23),  express  the  angles  in  radians,  then  to  obtain  au  in  the 
proper  quadrant,  rewrite  Equ.  (23)  in  the  form 


CK12  =  qatn 


A3 


/l  +  cain^y'*  /1  +  esin^y^ 

yl-csinfoy  ^1-csin^iy 


(24) 
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The  distance  is 


a« 


fa  —  ^1 


if  cosais  7^0, 


"cos  an'  (25) 

a«|  Aa  —  Ai| cos^i,    otherwise. 
In  the  WGS  1972  earth  model,  <z«  =  3443.9174  nautical  miles.  Also,  cr<ji  =  ai3  +  «r. 

b.  Direct  Solution.  Input  variables  are  the  latitude  fa  and  the  longitude  Ai  of  Pi , 
the  forward  aiimuth  aria  from  Pi  to  Pa,  and  the  distance  d  from  Pi  to  Pj.  The  output 
variables  are  the  latitude  fa  and  longitude  A?  of  Pa,  and  the  backward  azimuth  aai  from 
Pa  to  Pi .  All  angular  input  and  output  is  in  radians.  In  the  direct  rhumb  line  solution  for 
the  spheroid  earth,  if  cos  an  =  0,  then 


A,  =  Ai  +  ^?-,  and 


COB^j 


h-^u 


(26a) 


otherwise,  if  cos  aria  ^  0, 


Aa  =  Ai  ■+*  tan  oria 


/l  +  eaiiifay*      /l  +  gsinfrV 
.\1  —  tamfa)  ll  —  cainfa  J 


,  and 


(266) 


<h  =  <h+(d/a9)coacti2. 
Also,  orai  =  <*13  +  *• 

2.  Spherical  Earth  Formulas.  The  spherical  earth  formulas  are  obtained  by  setting 
t  -  0  in  Equations  (24)  and  (26b),  and  replacing  a«  by  the  standard  approximation  of  60 
nautical  miles  per  degree  of  arc  on  the  earth's  surface. 

a.  Inverse  Solution.  Input  variables  are  the  latitude  fa  and  longitude  Ai  of  Pi ,  the 
latitude  fa  and  longitude  Aa  and  P3,  fa  and  Aa.  Output  variables  are  distance,  d,  forward 
azimuth  from  Pi  to  Pa,  aia,  and  back  azimuth  from  Pa  to  Pi,  aai.  All  angular  input  and 
output  is  in  radians. 

The  spherical  earth  rhumb  line  course  is  obtained  by  setting  e  =  0  in  Equ.  (24).  The 
result  is 

<x„  =  qatn  L  -  Ai.ktan  ( J  +  fe)  -  htm  (|  +  it)! 

For  the  spherical  earth  model,  the  distance  d  in  nautical  miles  is 


(27) 


I      \*  )  cosoria' 
60/  —  ]  |Aj  -  Ai|co»^t,    otherwise. 
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(28) 


Also,  arai  =  aia  +  it. 

b.  Direct  Solution.  Input  variables  are  the  latitude  fa  and  the  longitude  Ai  of  Pi , 
the  forward  asimuth  a13  from  Pi  to  Pa,  and  the  distance  d  from  Pi  to  P?.  The  output 
variables  are  the  latitude  fa  and  longitude  A3  of  P3,  and  the  backward  asimuth  a^  from 
Py  to  Pi.  All  angular  input  and  output  is  in  radians. 
If  cosais  =0,  then 

*-Al  +  (lS>)(«)=?r,"d  (29a) 

otherwise,  if  cos  oris  ^  0, 

A3  =  Ai  +tanai3  [in tan  (f  +  y)  -In tan  (j  +  y)j>  and 

/w\/d\  (296) 

*3  =  *1  +  (li0J(60JCO8ai2- 

Also,  orji  =  cria  +  f. 
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m.  NAVALGOR:  Navigation  Algorithm  Program. 

A.  Introduction.  NAVALGOR  implements  three  sets  of  computational  procedures:  (1) 
Direct  Solution  Algorithms,  (2)  Inverse  Solution  Algorithms,  and  (3)  Rhumb  Line  Ap- 
proximations to  Great  Circle  Routes. 

The  "Direct  Solution  Algorithms*  compute  the  latitude  and  longitude  of  a  position  P3 
and  the  backward  azimuth  from  Pi  to  Pi  given  the  latitude  and  longitude  of  a  position  Pi , 
the  forward  asimuth  from  Pi  to  Pa  and  the  distance  from  Pi  to  Pj.  The  "Inverse  Solution 
Algorithms*  compute  the  distance  from  position  Pi  to  position  Pa,  the  forward  asimuth 
from  Pi  to  Pj,  and  the  backward  asimuth  from  Pa  to  Pi  given  the  latitude  and  longitude 
of  positions  Pi  and  P3.  For  comparision,  the  direct  and  inverse  computations  are  made 
using  four  procedures:  (1)  the  spheroid  earth  model,  (2)  the  spherical  earth  model,  and 
rhumb  line  approximations  using  both  (3)  the  spheroid  earth  and  (4)  the  spherical  earth 
models. 

The  spheroid  earth  model  should  be  considered  the  standard  of  comparison  for  all 
other  models.  Computations  performed  using  an  IBM  3033  and  double  precision  FORTRAN 
to  compare  the  inverse  solution  algorithm  [Equs.  (12)]  to  the  AGIO  long  lines  (50  to  6000 
miles)  [Ref.  9)  demonstrated  that  the  m«rimiiTfi  average  error  of  the  inverse  solution  algo- 
rithm was  -0.005  meters  with  a  standard  deviation  of  0.014  meters  (unpublished  result). 

The  spherical  earth  model  is  included  because  of  its  popularity  and  simplicity.  For 
many  procedures  the  spherical  earth  model  is  adequate,  and  it  certainly  is  more  compact 
for  use  in  a  small  computer  than  the  spheroid  earth  model. 

The  "Rhumb  Line  Approximations  to  Great  Circle  Routes9  procedure  may  be  used  to 
find  piecewise  constant  course  (rhumb  line)  approximations  to  great  circle  routes  for  any 
given  increment  in  longitude.  In  addition,  if  the  vertex  of  a  course  is,  for  example,  too  far 
to  the  north,  a  limiting  latitude  may  be  input  to  restrict  the  rhumb  line  approximation  to 
go  no  farther  north  than  that  limiting  latitude.  The  spheroid  earth  rhumb  line  equations 
are  used  in  this  section  of  the  program. 

In  the  sample  problems  and  in  the  program  listing,  the  east  minus  and  south  minus 
convention  is  used. 
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B.  Sample  Problems — NAVALGOR*  The  output  shown  in  these  sample  problems 
was  computed  using  Microsoft  basica/d  (double  precision)  on  an  IBM  PC. 
Master  Menu.  The  master  menu  for  NAVALGOR  is  shown  below. 

NAVIGATION  ALGORITHM  DEMO 

1)  DIRECT  SOLUTION 

2)  INVERSE  SOLUTION 

3)  RHUMB  LINE  APPROXIMATIONS 

TO  GREAT  CIRCLE  ROUTES 

4)  QUIT 

Option  (1)  selects  the  "Direct  Solution  Algorithm",  option  (2)  selects  the  "Inverse  Solution 
Algorithm",  option  (3)  selects  the  "Rhumb  Line  Approximations  to  Great  Circle  Routes" 
algorithm,  and  option  (4)  returns  the  user  to  the  operating  system. 

Problem  1.  Suppose  you  are  at  San  Francisco  (latitude  37°  47'  north  and  longitude 
122°25'  west),  that  your  initial  course  is  260°  and  that  you  travel  a  distance  of  4000  n.  mi. 
What  is  your  final  position?  Select  Option  1  from  the  master  menu. 

DIRECT  SOLUTION 

1st  LATITUDE      DD.MMSS  (-S)  ?  37.47 

1st  LONGITUDE  DDD.MMSS  (-E)  7  122.25 

INITIAL  COURSE    DDD.MMSS  7  260 

DISTANCE  (n.  ml.)  7  4000 

SPHEROID  EARTH  DIRECT  SOLUTION 
2nd  LATITUDE  6°  38.  7' 

2nd  LONGITUDE    -172°08.8' 
BACK  AZIMUTH         51°40.7' 

SPHERICAL  EARTH  DIRECT  SOLUTION 
2nd  LATITUDE  6°41.9' 

2nd  LONGITTJDE    -172°00.7' 
BACK  AZIMUTH         51°35.9' 

RHUMB  LINE  SOLUTIONS 

SPHERICAL  26°12.4/       -159°56.0' 

SPHEROID  28°12.4'      -160°18.4' 

PRESS  ANY  KEY  TO  CONTINUE 
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Problem  2.  Suppose  you  are  at  San  Francisco  (latitude  37°  47*  north  and  longitude 
122°25'  weat)  and  that  your  destination  is  Sydney,  Australia  (latitude  33°51'  south  and 
longitude  151°  13*  east).  How  far  do  you  travel,  what  is  your  initial  course,  and  what  is  the 
backward  asimuth  from  Sydney  to  San  Fransisco?  Select  Option  2  from  the  master  menu. 

INVERSE  SOLUTION 

1st  LATITUDE      DD.MMSS  (-S)  ?  37.47 

1st  LONGITUDE  DDD.MMSS  (-E)  7  122.25 

2nd  LATITUDE      DD.MMSS  (-S)  7  -33.51 

2nd  LONGITUDE  DDD.MMSS  (-E)  ?  -151.13 

SPHEROID  EARTH  INVERSE  SOLUTION 
DISTANCE  6443.52  n.al. 

FORWARD  COURSE        240°29.3' 
BACK  COURSE  55°55.7' 

SPHERICAL  INVERSE  SOLUTION 

DISTANCE  6446.3  n.al. 

FORWARD  COURSE        240° 18. 9' 
BACK  COURSE  55° 45. 9' 

RHUMB  LINE  SOLUTIONS:  SPHERE  SPHEROID 

DISTANCE  6464.49  6485.71  n.al. 

COURSE  228°19.7'       228°29.7' 

PRESS  ANY  KET  TO  CONTINUE 
Problem  3.  Suppose  you  are  at  latitude  37°  north,  longitude  76°  west  and  that  your 
destination  is  latitude  45°  north,  longitude  1°  west.  Compute  the  initial  great  circle  course 
and  distance,  the  latitude  and  longitude  of  the  vertex,  and  a  rhumb  line  approximation  to 
the  course  traveling  at  most  7°  degrees  of  longitude  on  each  leg  of  the  course.  In  addition, 
no  portion  of  the  course  is  to  be  more  northerly  than  47°.  Select  Option  3  from  the  master 
menu. 
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1st  Screens 

RHUMB  LINE  APPROXIMATIONS 

INITIAL  LATITUDE      DD.MMSS  (-S)  7  37 

INITIAL  LONGITUDE  DDD.MMSS  (-E)  7  76 

FINAL  LATITUDE  DD.MMSS  (-S)  ?  45 

FINAL  LONGITUDE      DDD.MMSS  (-E)  ?  1 

GREAT  CIRCLE  SOLUTION 

INITIAL  COURSE      56° 21. 2' 
DISTANCE  3307.8  n.mi. 

RHUMB  LINE  (MERCATOR)  SOLUTION 
COURSE  81° 58. 2' 

DISTANCE  3435.9  n.mi. 

VERTEX:  LATITUDE        48°19.8' 
LONGITUDE      28°07.3' 

PRESS  ANT  KEY  TO  CONTINUE 

In  this  first  screen,  the  great  circle  solution  and  a  single  rhumb  line  solution  are  computed 
and  displayed.  In  addition,  the  location  of  the  nearest  vertex  is  displayed.  The  vertex 
may  or  may  not  lie  between  the  initial  and  final  positions.  If  the  vertex  is  not  between 
the  initial  and  final  positions  then  any  limiting  latitude  which  may  be  input  on  the  next 
screen  prompt  will  be  ignored.  If  the  vertex  is  between  the  initial  and  final  positions,  a 
limiting  latitude  will  be  ignored  unless  it  is  lies  between  the  latitude  of  the  vertex  and  the 
latitude  of  the  position  which  is  closest  to  the  latitude  of  the  vertex. 
2nd  Screens 

RHUMB  LINE  APPROXIMATIONS 

INPUT  THE  MAXIMUM  NUMBER  OF  DEGREES 
OF  LONGITUDE  ON  EACH  RHUMB  LINE  LEG  ?  7 

LIMITING  LATITUDE  DD.MMSS  (-8) 

(0  TO  OMIT)  7  47 
LIMITING  LONGITUDES: 

10°  46. 8*      46°28.8' 

PRESS  ANY  KEY  TO  CONTINUE 

In  the  second  screen  a  prompt  appeared  for  the  maximum  longitude  between  each  rhumb 
line  leg — a  value  of  7°  was  input.   Also  input  was  a  limiting  latitude  of  47°,  which  lies 
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between  the  vertex  latitude  of  48°  19.8'  and  46° ,  the  latitude  of  the  final  position.  Limiting 
longitudes  of  10° 45.8'  and  45°28*8'  are  computed  and  displayed.  The  latitude  of  a  great 
circle  route  between  these  limiting  longitudes  will  be  more  northerly  than  47° ,  consequently 
any  rhumb  line  approximation  will  be  due  east  or  due  west  at  a  latitude  of  47°  between 
these  longitudes. 
3rd  Screenx 

RHUMB  LINE  APPROXIMATION  TO  GREAT  CIRCLE  COURSE 


GREAT 

GREAT 

RHUMB 

RHUMB 

CIRCLE 

CIRCLE 

LINE 

LINE 

LAT 

LONG 

COURSE 

DIST 

COURSE 

DIST 

37°  00.0* 

76°00.0/ 

56 

0 

58 

0 

39°27.9/ 

71°00.0' 

59 

278 

62 

279 

42°  18. 8' 

64°00.0' 

64 

639 

66 

641 

44°32.0/ 

57°00.0' 

69 

971 

71 

974 

46°  11. 7' 

50°00.0f 

74 

1283 

76 

1287 

47°00.0' 

45°28.8' 

77 

1475 

90 

1480 

47°00.0' 

10°45.8' 

103 

2884 

105 

2900 

45°56.3' 

S^O.C 

107 

3130 

108 

3148 

45°00.0' 

roo.c 

110 

3308 

108 

3326 

1)  NEW  APPROXIMATION  OR  2)  NEW  PROBLEM? 

Screen  3  displays  the  rhumb  line  approximation  for  a  maximnm  of  7°  of  longitude  between 
course  changes.  The  latitude  and  longitude  of  the  initial  and  final  positions  and  of  the 
positions  at  which  course  changes  are  made  are  given  in  columns  1  and  2,  respectively. 
Column  3  shows  the  great  circle  heading  at  each  position,  while  Column  4  shows  the 
cumulative  great  circle  distance  from  the  initial  position.  Similarly,  Column  5  shows  the 
rhumb  line  course  to  be  followed  between  each  pair  of  positions  and  Column  6  shows 
the  cumulative  rhumb  line  distance  from  the  initial  position.  Note  that  at  47°  north  and 
46°  28.8'  west  (the  first  limiting  latitude  position)  the  rhumb  line  course  changes  to  due  east 
until  47°  north  and  10°  45.8'  west  in  reached.  Also  note  that  the  rhumb  line  approximation, 
including  the  'detour1  at  the  limiting  latitude,  is  only  18  n.  mi.  longer  than  the  great  circle 
route.  The  third  screen  ends  with  two  options:  Option  1  returns  to  Screen  2  and  prompts 
for  new  inputs  while  Option  2  returns  to  the  master  menu. 
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C.  Program  Listing-- NAVALGOR. 

10  BEN  "NAVIGATION  ALGORITHMS"  R.H.  SHUDDE,  03-12-85.  REV.  00-14-85  1600 

13  REN  "NAVALGOR/BAS- 

30  DEFDBL  A-Y 

40  P4-ATNC1) :P>P4*P4:PI-P2*P2:TP-PI+PI:RD"PI/180:EP-lE-33 

50  FW/2©8.26:AE-637813S/1862:E<>SQR(FL*(2-FL)) 

60  DEFFNN(X)-X-MO*INT(X/MO):REMX  MOD  MO  FUNCTION. 

70  DEFFNL(X)-X-TP*INT((X*PI)/TP):REN  LONGITUDE  ADJUST  (-PI.PI) 

80  DEFFNR(X)-INT(X*N0+.5)/M0:REM  ROUNDING  FUNCTION. 

00  DEFFNG(X)-X+PI*8GN(X)*(AB8(X)>P2):REM  UTITUDE  ADJUST  (-PI/2.PI/2) 

100  DEFFNC(I)-ATN(SqR(l-X*X)/(X-EP*(X-0#)))-PI*a«'0#):REM  ARCC08 

110  DEFFNSa)=ATN(X/(SqR(l-X*X)-EP*(ABS(X)=l))):REM  ARCSIN 

115  GOTO  2000 

120  REN  QATN  (-PI.PI)  FUNCTION: 

130  A»ATN(Y/(X-EP*(X=0#)))-PI*(X<0#)*(SGN(T)-(Y=Oi)) :RETURN 

140  REN  QATN  (O.TWOPI)  FUNCTION: 

150  A-ATNCY/  (X-EP*  (X=0#)  )  )  -PI*  (X<0#)  +TP*  (X>-0#)  *  (Y<0#)  :  RETURN 

170  REN 

200  REN  DIRECT  SOLUTION,  SPHEROID  EARTH.  ALL  ANGLES  MUST  BE  IN  RADIANS. 

210  REM  INPUT:  UTITUDE  Gl,  LONGITUDE  LI.  FORWARD  AZIMUTH  Al  AND 

220  REN  DISTANCE  DD*(S/AE)  TO  A  POINT  P2.  NOTE:  DD  HAS  RADIAN  UNITS. 

230  REN  OUTPUT:  UTITUDE  G2.  LONGITUDE  L2  AND  BACKWARD  AZIMUTH  A2. 

240  S0-8IN(A1):C0-C08(A1):TA-ATN((1-FL)*TAN(G1)) 

250  88-SIN(TA):C8-C0S(TA):N—  SQ*C8:C1«FL*N 

260  OFL*(l-M*M)/4:D-a-C2)*(l-C2-Ci*M) 

270  P-C2*(1+.5*C1*N)/D:N-C8*C9 

280  Y-N:X-S8:G0SUB  130:8G*A:8D-DD/D:U*2*(SG-8D) 

200  W-1-2*P*C0S(U) :V=COS(U+SD) :X=C2*C2*8IN(SD)*C08(SD*(2*V*V-1)) 

300  Y-2*P*V*W*SIN(SD)  :D0-8D*X-Y:SL-8IN(Dq)  :CL-COS(DQ) 

310  X-8qR(M*M+(N*CL-88*8L)A2)  :G>ATN((S8*CL+N*SL)/X/(1-FL)) 

320  Y— 8L*S9:X«C8*CL-S8*SL*C8:G0SUB  130:DJ*A 

330  B-Cl*(l-C2)*Dq-Cl*C2*8L*C0S(SG*SG-0q) 

340  L2=FNL(L1+DJ-H):Y»N:X«-(N*CL-88*SL):G0SUB  150:A2*A 

360  RETURN 

360  REM 

400  REM  INVERSE  SOLUTION,  SPHEROID  EARTH.  ALL  ANGLES  MUST  BE  IN  RADIANS. 

410  REM  INPUT:  UTITUDES  Gl  *  G2,  AND  LONGITUDES  LI  *  L2. 

420  REM  OUTPUT:  DISTANCE  DD»(S/AE).  FORWARD  AZIMUTH  Al,  AND  BACKWARD 

430  REM  AZIMUTH  A2.  NOTE:  DD  HAS  RADIAN  UNITS. 

440  TA-ATN(C1-FL)*TAN(G1)):TB«ATN((1-FL)*TAN(G2)) 

450  DN-L2-Ll:DT-(TB-TA)/2:TM-(TA*TB)/2 

460  SH-8IN(DT):CH-C08(DT):8M-8IN(TM):CM-C08(TM) 

470  H-CH*CH-8M*SM:L-8H*SH+H*(SIN(DM/2))~2 

480  SD-2*FNS(SqR(L)):U-2*8M*SM*CH*CH/((l-L)-(EP*(L-l))) 

400  V«2*8H*SH*CM*CM/L:X-U>V: Y-U-V: T-8D/SIN(SD) 

500  D-4*T*T:E-2*C0S(SD) :A-D*E:C-T-(A-E)/2:N1-X*(A-K>X) 

510  B-D+D : N2-Y* (B+E*Y) : N3-D*X*Y : D4-FL*FL* (N1-N2+N3) /64 

520  D3-FL*(T*X-Y)/4:DD-(T-D3+D4)*SIN(SD) 

530  M-32*T- (20*T-A) *X- (B+4) *Y : F=Y+Y-E* (4-X) 

540  G-FL*(T/2+FL*M/64) :Q— (F*G*TAN(DM))/4 

550  DW-(DM+q)/2:SW-8IN(DW):CW-C08(DW) 
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560  Y*8H*CW:X-CM*8W:G08UB  130:T1-A 

670  Y— CH*CI:X»8M*8¥:G0SUB  130:T2=A 

680  MO^P:Al-FNM(Tl+T2):A>FNM(Tl-T2):RETUHH 

600  REM 

000  REM  DIBECT  SOLUTION,  SPHERICAL  EABTH.  ALL  ANGLES  MUST  BE  IN  RADIANS. 

010  REN  INPUT:  LATITUDE  61,  LONGITUDE  LI.  FORWARD  AZIMUTH  Al  AND 

620  REN  DISTANCE  DD  TO  A  POINT  P2.  NOTE:  DD  IS  IN  RADIANS. 

630  REM  OUTPUT:  LATITUDE  G2.  LONGITUDE  L2  AND  BACKWARD  AZIMUTH  A2. 

640  S9=8IN(Al):CO=C0S(Al):S3=8IN(Gl):C3=C0S(Gi) 

660  SD-SIN(DD): CD-COS (DD) 

660  Y«8D*SO:X<3*CD-B3*C0*SD:G0SUB  130:L2*FNL(L1-A) 

670  Y— S©*C3:X=»8D*83-CD*C3*C3;G0SUB  1S0:A2»A 

680  G2-FNS(S3*CD+C3*SD*C0) :  RETURN 

600  REM 

700  REN  INVERSE  SOLUTION,  SPHERICAL  EARTH.  ALL  ANGLES  MUST  BE  IN  RADIANS. 

710  REM  INPUT:  LATITUDES  Gl  ft  G2.  AND  LONGITUDES  LI  ft  L2. 

720  REM  OUTPUT:  DISTANCE  DD  TO  A  POINT  P2.  (NOTE:  0  <=  DD  <=»  PI  RADIANS) 

730  REM  FORWARD  AZIMUTH  Al.  AND  BACKWARD  AZIMUTH  A2. 

740  S3-8IN(G1) :C3-C0S(G1) :S4-8IN(G2) :C4«C08(G2) 

760  DL»L1-L2:SU-8IN(DL):CU=C08(DL) 

760  DD-FNC(S3*84+C3*C4*CU) 

770  Y-8U*C4:X-C3*84-83*C4*CU:G08UB  150:A1-A 

780  Y— 8U*C3:X»C4*S3-S4*C3*CU:G08UB  160 :A2=A: RETURN 

700  REM 

800  REM  INVERSE  SOLUTION.  RHUMB  LINE  FOR  SPHERICAL  ft  SPHEROIDAL  EARTH. 

810  REM  INPUT:  LATITUDES  Gl  ft  G2.  AND  LONGITUDES  LI  ft  L2. 

820  REM  OUTPUT:  DISTANCES  DA  ft  DB,  AND  AZIMUTHS  ZA  ft  ZB. 

830  GA=TAN(P4+Gl/2) :GB=TAN(P4+G2/2) :DL=FNL(L1-L2) : DK=ABS(DL/RD) 

840  GA-GA-EP* (GA-O) : GB-GB-EP* (GB-O) 

860  DG»(G2-Gl)/RD:Y=DL:X»L0G(GB)-L0G(GA):GOSUB  150:2A=A 

860  CZ-C08(ZA):DA-60*DE*COS(G1):IF  AB8(CZ)>. 0001 THEN  DA-60*DG/CZ 

870  El-EC*SIN(Gl):E2-EC*SIN(G2):ED-EC/2 

880  Y-DL:X-L0G(GB/((1*E2)/(1-E2))-ED)-L0G(GA/((1+E1)/(1-E1))*ED) 

886  GOSUB  160:ZB-A 

800  CZ-COS<ZB):DB-60*DI*CQS(G1):IF  AB8(CZ)>  0001 THEN  DB-60*DG/CZ 

000  RETURN 

010  REN 

1000  REN  DIRECT  SOLUTION.  RHUMB  LINE  FOR  SPHERICAL  ft  SPHEROIDAL  EARTH. 

1010  REN  INPUT:  LATITUDE  Gl.  LONGITUDE  LI.  FORWARD  AZIMUTH  Al  AND 

1020  REN  DISTANCE  DD  TO  A  POINT  P2.  NOTE:  DD  IS  IN  RADIANS. 

1030  REN  OUTPUT:  LATITUDE  G2  AND  LONGITUDE  LR  ft  LS 

1040  C6-CQ8(A1) 

1060  IF  C8=0THEN  LR=Ll+DD/C08(Gl):LS=La:G2=Cl: RETURN 

1060  G2-G1+DD*C8 

1070  GA-TAN(P4+Gl/2) :GB-TAN(P4+G2/2) : T6-TAN(A1) 

1080  LR-L1-T6*(L0G(GB)-L0G(GA)) 

1000  El-EC*8IN(Gl):E>EC*8IN(G2):ED-EC/2 

1100  LS»Ll-Te*(L0G(GB/((l4E2)/(l-E2))-ED)-LOG(GA/((l+El)/(l-El))'*ED)) 

1110  RETURN 

1120  REN 

1200  REN  DECIMAL  TO  DDD  MN.F 

1210  V$»«  ":IF  KOTHEN  V$-«--:X— X 

1220  X«X-H/1200:Y-INT(X):V$-V$*RIGHT$("  ■+8TR$(Y),3)*''0" 
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1230  X-600*(X-Y):T-INT(X):X|-STR|<1000*Y) 

1240  V|-V|*MIDiail3,2)*-.«+RICHT$(X|ll)+"'":aErnBN 

1260  REM 

1280  BEN  DDD.MMS8  TO  DECIMAL 

1270  IX-0:P01  >1T0  LZMCVDjCI-MIDKVI.Z.DiIF  C|-"."THEN  IX-Z 

1280  NEXT  :IF  IX-OTHEN  X-VAL(Vl)  : RETURN 

1200  X-VALCLEFTKVI.ID^SN*!:!?  KOTHEN  8M— 8N:X— X 

1300  Vt-V$^0000«:T*VAL(MID$(Y$>n*1.2)):Z-VAL(MIDt(Vt.IX*3,2)) 

1310  X-8N*((Z/60*Y)/60+X):RETURN 

1320  REM 

2000  CLS  : PRINT  8PC( 15) ; "NAVIGATION  ALGORITHM  DEMO 

2010  PRINT  : PRINT  : PRINT 

2020  PRINT  SPC<15);"1)  DIRECT  SOLUTION 

2030  PRINT  SPC(1S);"2)  INVERSE  SOLUTION 

2040  PRINT  SPC(1S);'3)  RHUMB  LINE  APPROXIMATIONS" 

2060  PRINT  SPC(iS) ; •  TO  GREAT  CIRCLE  ROUTES 

2060  PRINT  :PRINT  SPC(15);"4)  QUIT 

2070  GOSUB  6010:C=VAL(C|):ON  CG08UB  3000,3600,4000,5500 

2080  GOTO  2000 

2000  REM 

3000  CLS  : PRINT  SPC(15); "DIRECT  SOLUTION" .PRINT  : PRINT 

3010  PRINT  8PC(10);: PRINT  "1ST  LATITUDE  DD.MM88  (-8)  "; 

3020  INPUT  VI: GOSUB  1270:G1-X*RD 

3030  PRINT  SPC(IO);: PRINT  "1ST  LONGITUDE  DDD.MM88  (-E)  "; 

3040  INPUT  VI: GOSUB  1270:L1-X*RD 

3060  PRINT  SPC(10);:PRINT  "INITIAL  COURSE  DDD.MM8S  "; 

3060  INPUT  VltGOSUB  1270:A1-X*RD 

3070  PRINT  SPC(IO); "DISTANCE  (N.MI.)  ";:INPUT  D 

3080  D1-D*RD/60:DD-D/AE 

3000  GOSUB  240: PRINT  : PRINT  SPC (8) J "SPHEROID  EARTH  DIRECT  SOLUTION 

3100  MO-100:  GOSUB  3210 

3110  DD-D1: GOSUB  640:PRINT  : PRINT  SPC (8) ; "SPHERICAL  EARTH  DIRECT  SOLUTION 

3120  GOSUB  3210 

3130  GOSUB  1040: PRINT  : PRINT  SPC(8) ; "RHUMB  LINE  SOLUTIONS"; 

3136  PRINT  8PC(4);"LAT";SPC(14);"L0NG 

3140  PRINT  SPC(12) ; "SPHERICAL  "; 

3150  X-FNG(G2):X»X/RD:G0SUB  1210 :VSI-V|: PRINT  V|; SPC (8); 

3160  X-FNL(LR) :X«X/RD: GOSUB  1210:PRINT  VI 

3170  PRINT  SPCC12) ; "SPHEROID  ";V8|;SPC(8); 

3180  X-FNL(LS):X-X/RD: GOSUB  1210 -.PRINT  V| 

3100  GOSUB  6000: GOTO  2000 

3200  REM 

3210  PRINT  SPC(12);"2ND  LATITUDE  "; :X"Q2/RD: GOSUB  1210:PRINT  V| 

3220  PRINT  8PC(12);"2ND  LONGITUDE  •; :X»L2/RD: GOSUB  1210:PRINT  VI 

3230  PRINT  SPCQ2);"BACX  AZIMUTH  ■;: X-A2/RD : GOSUB  1210:PRINT  V|: RETURN 

3240  REM 

3600  CLS  :PRINT  8PC(15); "INVERSE  SOLUTION" : PRINT  :PRINT 

3610  PRINT  SPC(10);:PRINT  "1ST  LATITUDE  DD.MMS8  (-8)  ■; 

3620  INPUT  VI: GOSUB  1270:G1-X 

3630  PRINT  8PC(10); -.PRINT  "1ST  LONGITUDE  DDD.MM88  (-E)  ■; 

3640  INPUT  VI: GOSUB  1270:L1-X 

3660  PRINT  SPC(IO);: PRINT  "2ND  UTITUDE  DD.MMSS  (-S)  ■; 

3660  INPUT  VI: GOSUB  1270 :G2-X 

24 


3570  PRINT  SPC(IO);: PRINT  "2ND  LONGITUDE  DDD.MMSS  (-E)  ■; 

3680  INPUT  V$:GOSUB  1270:L2-X 

3600  G1-G1*RD:G2-G2*RD:L1-L1*RD:L2-L2*RD 

3600  OOSUB  440: PRINT  : PRINT  SPC(8) ; •SPHEROID  EARTH  INVERSE  SOLUTION 

3610  MO- 100:11-1: GOSUB  3700 

3620  DD-D1:G0SUB  740:PRINT  :PRINT  SPC (8) ; "SPHERICAL  INVERSE  SOLUTION 

3630  IX=2:G0SUB  3700:G08UB  830 

3640  PRINT  : PRINT  SPC (8) ; -RHUMB  LINE  SOLUTIONS:  SPHERE  SPHEROID 

3660  PRINT  8PCC12); "DISTANCE  ";FNR(DA);3PC(8);FNR(DB);"  N.MI. 

3660  PRINT  SPC(12); "COURSE  "; :X-ZA/RD:60SUB  1210:PRINT  V|;SPC(6); 

3670  X=ZB/RD:GOSUB  1210:PRINT  V| 

3680  GOSUB  6000: GOTO  2000 

3600  RSI 

3700  PRINT  SPCC12); "DISTANCE  "; 

3710  IF  IX'ITHEN  PRINT  FNR(DD*AE);"  N.MI. 

3720  IF  IX-2THEN  PRINT  FNR(60*DD/RD) ; ■  N.MI. 

3730  PRINT  SPC(12);"F0R»ARD  COURSE  ■; :X-A1/RD: GOSUB  1210:PRINT  V| 

3740  PRINT  SPC(12);"BACX  COURSE  ■; :X-A2/RD: GOSUB  1210:PRINT  V|:RETURN 

3760  REM 

4000  CLS  :  PRINT  SPC(  16)  ;  "RHUMB  LINE  APPROXIMATIONS" :  PRINT  :  PRINT 


4010  PRINT  SPC(IO); 

4020  INPUT  VI :GOSUB 

4030  PRINT  SPC(IO); 

4040  INPUT  V$:GOSUB 

4060  PRINT  SPC(IO); 

4060  INPUT  V$:GOSUB 

4070  PRINT  SPC(IO); 


; PRINT  "INITIAL  LATITUDE  DD.MM83  (-8)  "; 

1270:G1-X*RD:06-01 

: PRINT  "INITIAL  LONGITUDE  DDD.MMS8  (-E)  "; 
1270:L1»X*RD:L6-L1 

:PRINT  "FINAL  LATITUDE  DD.MMS8  (-8)  "; 
1270:G2-X*RD:G7-G2 

: PRINT  "FINAL  LONGITUDE  DDD.MMSS  (-E)  "; 
4080  INPUT  7$:  GOSUB  1270 :  L>X*RD :  L7-L2 

4000  GOSUB  740 :M0=10: PRINT  : PRINT  SPC ( 10) ; "GREAT  CIRCLE  SOLUTION 
4100  PRINT  SPC(IS); "INITIAL  COURSE  ■; :AZ-A1:X-A1/RD : GOSUB  1210:PRINT  71 
4110  PRINT  8PCC15); "DISTANCE  •;FNR(60*DD/RD);"  N.MI. 
4120  GOSUB  330:PRINT  : PRINT  SPC(12); "RHUMB  LINE  (MERCATOR)  SOLUTION 
4130  PRINT  SPC(IS); "COURSE  ■; :X»ZB/RD: GOSUB  1210:PRINT  VI 
4140  PRINT  SPC(IS); "DISTANCE  ";FNR(DB);«  N.MI. 
4160  GOSUB  4610 

4160  PRINT  :PRINT  SPC(12); "VERTEX:  LATITUDE  "; :X-GV/RD: GOSUB  1210:PRINT  V$ 
4170  PRINT  8PC(20) ; "LONGITUDE  ■; :X-LV/RD: GOSUB  1210: PRINT  V$ 
4180  GOSUB  6000 

4190  CLS  :PRINT  SPC(16); "RHUMB  LINE  APPROXIMATIONS": PRINT 
4200  PRINT  :  PRINT  SPC(IO);  "INPUT  THE  MAXIMUM  NUMBER  OF  DEGREES 
4210  PRINT  SPC(IO);"  OF  LONGITUDE  ON  EACH  RHUMB  LINE  LEG  ■; 
4220  INPUT  VI: GOSUB  1270:IN-X*RD 
4230  GS-0:R8-0:IF  IN-OTHEN  IN-TP 
4240  DL-FNLCL7-L6) : IF  DL-OTHEN  F-l 
4260  DV-FNLCLV-L6)  :  IF  DLoOTHEN  F-DV/DL 
4260  LX-0:IF  F<-O0R  F>-1THEN  LX-1:G0T0  4370 
4270  PRINT  : PRINT  SPC(IO); "LIMITING  UTITUDE  DD.MM8S  (-8) 
4280  PRINT  SPC(IO);"  (0  TO  OMIT)  ";:INPUT  V|:G08UB  1270:LL-X*RD 
4290  L8-0 :  L9-0 :  G1-G8 :  G3-G1 :  L1-L6 :  Ul-Ll :  PH-LL :  G2-G7 :  L2-L7 
4300  IF  LL-OTHEN  LX-1:G0T0  4370 

4310  Q1-AB8(FV-G6)  :  q>AB8(PV-G7)  :  A-qi :  IF  Q2<A  THEN  A-Q2 
4320  A-A*8GN(PV):B-PV-LL:IF  A-0  THEN  G-6:G0T0  4340 
4330  G-B/A 
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4340  IF  6<-0  OR  G>1  TIEN  I>1:P&INT  SPC(10);"CANT  USE.  -  IGNORED. ■: GOTO  4370 

4360  GOSUB  4810: PRINT  SPC(IO); 'LIMITING  LONGITUDES: 

4360  X-L8/RD:G0SUB  1210:PRINT  SPCC13);V$;SPC(10):>L9/RD:G0SUB  1210:PRINT  V$ 

4370  GOSUB  6000: GOSUB  5220 

4380  D8-FNL(L8-L1):D>PNL(L0-L1):IP  AB8(D0)<ABS(D8)THEN  T-L9:L9-L8:L8-T 

4300  MX-1:LZ-L6:PX-G6:IF  LX-1THEN  4460 

4400  LT-LB:PY«LL: GOSUB  4010 

4410  Li«L8:Gl»LL:L2»L9:G2«LL: GOSUB  740:G0SUB  830 

4420  GOSUB  5110:GS-G8+DD:RS-R8+DB 

4430  LX-L9:PX-LL:L1-LX:01-PX 

4440  L2»L7  G2=G7 

4460  MX-0:LY-L7:PY-G7: GOSUB  4010 

4460  REM 

4470  PRINT  : PRINT  SPC(10);N1)  NEW  APPROXIMATION  OR  2)  NE*  PROBLEM? 

4480  GOSUB  6010: OVAL(C$)  : ON  CGOSUB  4100,2000 

4400  GOTO  4480 

4500  REM 

4600  REM  COMPUTE  VERTEX 

4610  8A-8INCAZ) :CA-C08(AZ) :SP-8IN(G1) :CP-C0S(G1) 

4620  IF  3A-0THDJ  GV»P2:LV-L8: RETURN 

4630  LY-L1-ATN(1/SP/TAN(AZ)):U-L1:LV-FNL(LV):LM-LV:G0SUB  4720 

4640  GY=PH:  RETURN 

4660  REM 

4700  REM  FIND  LATITUDE  GIVEN  LONGITUDE 

4710  SA-8IN(AZ) :CA-COS(AZ) :SP-8IN(PA) :CP-COS(PA) 

4720  DL-LA-LM:Y-SP*C0S(DL)*SA+8IN(DL)*CA:X-CP*SA:  GOSUB  130:PI-A 

4730  PH-FNG(PI): RETURN 

4740  RSI 

4800  REN  FIND  LONGITUDE  GIVEN  LATITUDE 

4810  KX-0:IF  ABS(PH)->AB8(GV)THEN  KX-1: RETURN 

4820  8A-SIN(AZ):CA-C08(AZ):SP-8IN(G3):CP-C0S(G3) 

4830  K-CP*8A*TAN(PH):(>8P*8A:S-CA:RQ-8qR(S*8-K:*C):Y-8:X-C:G0SUB  130:NU-A 

4840  CC-PNC(X/RO)  :DD-LA1-NU:L8-FNL(DD-CC)  :L9-FNL(DD*CC)  : RETURN 

4850  REM 

4900  REM  RHUMB  APPROXIMATIONS  SUBROUTINE.  INPUT:  LX,  PX.  LY  *  PT. 

4010  DL=FNL(LY-LX)  :  IOABS(IN)  :  AG=ABS(DL)/IC 

4920  IF  IC/RD<1THEN  NP-INTCAG>.5):M0».l 

4930  IF  IC/RD>*1THEN  NP*INT(AG) :M0*1 

4940  IF  NP<1TBEN  IC=DL : L2=LY i GOTO  4970 

4960  IC-SGN(DL)*IC 

4960  L2-FNR((LX+(DL-(NP-l)*IC)/2)/RD+.5*SGN(DL))*RD 

4970  L1=LX:G1=PX;LM=L2:GGSUB  4720:G2»PH 

4980  GOSUB  740G08UB  830GOBUB  5110:GS=CS*DD:RS=RS+DB 

4990  IF  NP<1THEN  5060 

5000  NOliIF  NP-1THEN  5040 

5010  L1-L2 : G1-G2 : L2-FNL (Li +IC) : LM-L2 : GOSUB  4720 : G2-PH 

5020  GOSUB  740: GOSUB  830: GOSUB  5110:GS=GS+DD:RS=RS+DB 

5030  N0N01:IF  NC<NPTHEN  5010 

5040  L1-L2:G1-G2:L2-LY:C2-PY 

5060  GOSUB  740:G0SUB  830:G08UB  5ilO:G8-GS+DD:RS-RS+DB:IF  MX-1THEN  RETURN 

5060  L1-LY:G1-PY:M0-TP:A1-FNM(A2*PI): GOSUB  51 10 -.RETURN 

5070  REM 

5100  REM  PRINT  LINE  OF  OUTPUT 
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5110  M0-1:X-G1/RD:00SUB  1210:PRINT  V$;SPC(2); 

5120  X-L1/BD:G0SUB  1210:PRINT  V$; 

5130  X-FNRCA1/RD) :  GOSUB  5200 :  PRINT  V$ ;  SPC(4)  J 

6140  X=FNR(00*GS/RD):GQSUB  5200:PRINT  V$;SPC(4); 

5150  X-FNRCZB/RD): GOSUB  5200:PRINT  Y$;SPC(3); 

5160  X»FNR(RS): GOSUB  6200:PEINT  V$ 

5170  NL~NL+1:IF  NL=20THEN  G08UB  6000: GOSUB  5220 

5180  RETURN 

6100  REM 

6200  V$-«  g+8TRia):V$-RIGHT$(V$, 7)  .RETURN 

5210  REN 

6220  CLS  : PRINT  SPC(IO) ;  "RHUMB  LINE  APPROXIMATION  TO  GREAT  CIRCLE  COURSE1 

PRINT 

5230  PRINT  ■  GREAT  GREAT  RHUMB  RHUMB 

6240  PRINT  '  CIRCLE  CIRCLE  LINE  LINE 

5250  PRINT  "  LAT  LONG  COURSE  DI3T  COURSE  DIST 

5260  NL-6: RETURN 

5270  REM 

6500  CLS  :END 

5510  REM 

6000  PRINT  : PRINT  SPC(10);aPRE88  ANY  KEY  TO  CONTINUE 

6010  FOR  I1-1T0  d:C*-INXEY$  :NEXT 

6020  d-DHEY*  :IP  C$— TREN  6020 

6030  RETURN 
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D.  Program  Annotation— NAVALGOR* 

Line(s)  Usage 

40  PI4  =  */4.  PI2  =  t/2.  PI  =  *.  TP  =  2*.  RD  =  */180  is  the  degree-to-radian 

conversion. 
50  FL  is  the  earth's  flattening  factor.  AE  is  the  earth's  equatorial  radius.  (Both 
of  these  constants  are  from  the  WGS-72  earth  model.)  ECC  is  the  earth's 
eccentricity. 
60  FNM  is  the  x  mod  MO  function. 
70  FNL  adjusts  longitude  to  lie  between  -f  and  w. 
80  FNR  rounds  x  to  the  nearest  M0* h. 
90  FNG  adjusts  latitudes  to  lie  between  -ir/2  and  ir/2. 
100  FNACS  is  the  arccosine  function.  In  the  Sharp  PC-1500A,  function  calls  may 

be  replaced  inline  by  the  ACS  function. 
110  FNASN  is  the  arcsine  function.  In  the  Sharp  PC- 1500 A,  function  calls  may  be 

replaced  inline  by  the  ASN  function. 
130  FNATN2  is  the  qatn  function  which  returns  a  principle  value  between  -*■  and 

w. 
140  FNATNP  is  a  qatn  function  which  returns  a  principle  value  between  0  and  2*°. 
160  GOTO  2000  to  commence  execution. 
240-350  Computation  of  the  spheroid  earth  direct  solution.  Gl  =  fa,  Ll  =  A1}  Al  = 
oeiq  and  distance  DD  are  input.  G2  =  fa,  L2  =  Ag  and  A2  =  a^\  are  output. 
440-580  Computation  of  the  spheroid  earth  inverse  solution.  Gl  =  fa ,  Ll  =  Ai ,  G2  = 
fa  and  L2  =  \i  are  input.  Al  =  atia,  A2  =  a?i  and  distance  DD  are  output. 
640-680  Computation  of  the  spherical  earth  direct  solution.  Gl  =  ^i,  Ll  =  Ai,  Al  = 
an  and  distance  DD  are  input.  G2  =  fa,  L2  =  Aq  and  A2  =  a?i  are  output. 
740-780  Computation  of  the  spherical  earth  inverse  solution.  Gl  =  ^i,  Ll  =  Ai,  G2  = 
fa  and  L2  =  A2  are  input.  Al  =  ar13,  A2  =  cr^i  and  distance  DD  are  output. 
830-900  Computation  of  the  rhumb  line  inverse  solution.  Gl  =  <^i,  Ll  =  Ai,  G2  =  fa 
and  L2  =  Aq  are  input.  Spherical  earth  distance  DA  and  azimuth  ZA  as  well 
as  spheroid  earth  distance  DB  and  azimuth  ZB  are  output. 
1040-1100  Computation  of  the  rhumb  line  direct  solution.    Gl  =  <^i,  Ll  =  A1,Al  = 
or  12  and  distance  DD  are  input.  G2  =  fa,  spherical  earth  longitude  LR  and 
spheroid  earth  longitude  LS  are  output. 
1200-1240  Convert  decimal  degrees  to  degrees,  minutes  and  tenths  of  minute  format. 
1270-1310  Convert  packed  degrees,  minutes  and  seconds  format  (DDD.MMSS)  to  decimal 

degrees. 
2000-2060  Master  menu  option  display. 

2070  Transfer  to  selected  menu  item. 
3000-3230  Direct  solution  option. 

3080  D  'is  the  distance  (input)  in  nautical  miles,  Dl  is  the  spherical  earth  distance 

in  radians  and  DD  is  the  distance  in  earth  equatorial  radius  units. 
3090  Compute  spheroid  earth  direct  solution.  Print  heading. 
3100  Print  spheroid  earth  direct  solution. 
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3110 

3120 

3130 

3140-3160 

3170-3180 

3190 

3210-3230 

3600-3740 

3600 

3610 

3620 

3630 

3640 

3650 

3660-3670 

3680 

3700-3740 

4000-5260 

4000-4080 

4090-4110 

4120-4140 

4150-4170 

4180 

4190-4220 

4230 


4240-4260 

4270-4280 
4290-4340 


4350-4360 


4370 
4380 

4390 


Compute  spherical  earth  direct  solution.  Print  heading. 

Print  spherical  earth  direct  solution. 

Compute  rhumb  line  solution.  Print  heading. 

Print  spherical  earth  rhumb  line  latitude  and  longitude. 

Print  spheroid  earth  rhumb  line  latitude  and  longitude. 

Continue  prompt.  Go  to  master  menu. 

Output  subroutine  for  option  1. 

Inverse  solution  option. 

Compute  spheroid  earth  inverse  solution.  Print  heading. 

Print  spheroid  earth  inverse  solution. 

Compute  spherical  earth  inverse  solution.  Print  heading. 

Print  spherical  earth  inverse  solution. 

Heading  for  rhumb  line  solutions. 

Print  spherical  and  spheroid  earth  rhumb  line  distances. 

Print  spherical  and  spheroid  earth  rhumb  line  courses. 

Continue  prompt.  Go  to  master  menu. 

Output  subroutine  for  option  2. 

Rhumb  line  approximation  to  great  circle  option. 

Input  prompts. 

Compute  great  circle  solution  and  print  initial  course  and  distance. 

Compute  rhumb  line  solution  and  print  course  and  distance. 

Compute  and  print  the  latitude  and  longitude  of  the  vertex  of  the  great  circle 

route. 

Display  continue  prompt. 

Prompt  for  the  longitude  increment  of  the  rhumb  line  approximation. 
GS  and  RS  are  the  cumulative  great  circle  and  rhumb  line  distances  traveled  on 

each  leg,  respectively.  Initialise  them  to  zero.  IN  is  the  longitude  increment 

in  radians.  If  input  as  zero  (no  increments  requested),  it  is  set  to  2*\ 
Determine  if  the  vertex  lies  between  the  origin  and  destination.  If  not,  go  to 

4370.  Otherwise  proceed. 

If  the  vertex  is  on  the  great  circle  route,  prompt  for  a  limiting  latitude  LL. 
Set  limiting  longitudes  L8  and  L9  equal  to  zero.    Determine  if  the  limiting 

latitude  cuts  the  great  circle  course.  If  it  does  not,  inform  the  user  that  the 

limiting  latitude  is  ignored.  L%  is  zero  if  the  limiting  latitude  is  to  be  used, 

otherwise  it  is  one. 
If  the  limiting  latitude  cuts  the  great  circle  course,  compute  and  display  the 

longitudes  L8  and  L9  at  which  the  limiting  latitude  cuts  the  great  circle. 
Display  continue  prompt,  then  print  heading. 
Determine  which  limiting  longitude  is  closest  to  the  initial  position.  L8  becomes 

the  closest,  L9  the  farthest. 
LX  and  PX  are  the  initial  great  circle  longitude  and  latitude.  If  the  limiting 

latitude  is  not  to  be  used,  go  to  4450. 
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4400  LY  and  PY  are  the  longitude  and  latitude  of  the  final  point  on  the  first  segment 
of  the  great  circle  course,  gosub  4910  to  compute  the  rhumb  line  approxima- 
tion up  to  the  first  limiting  longitude. 
4410  Compute  the  great  circle  and  rhumb  line  distances  from  the  limiting  latitude 

at  the  first  limiting  longitude  to  the  second  limiting  longitude. 
4420  Print  one  line  of  output  (for  the  limiting  longitudes).    Update  the  distance 
counters. 
4430-4440  Reset  LX,  PX,  LY  and  PY  for  the  section  of  the  great  circle  and  rhumb  line 
following  the  limiting  latitude  section  of  the  course. 
4450  GOSUB  4910  to  compute  the  final  section  of  a  course  following  a  limiting  latitude 
leg  or  compute  the  entire  course  for  cases  in  which  the  limiting  latitude  is  not 
used. 
4470-4490  Prompt  for  either  a  rework  of  the  rhumb  line  approximation  with  new  param- 
eters or  for  a  new  problem. 
4610-4640  Compute  the  latitude  GV  an  longitude  LV  of  a  great  circle  vertex. 
4710-4730  For  a  given  longitude,  compute  the  corresponding  latitude  of  on  a  great  circle. 
4810-4840  For  a  given  latitude,  compute  the  corresponding  longitudes  L8  and  L9  on  a 
given  great  circle.  K%  is  set  to  one  if  there  is  no  solution,  otherwise  K%  is 
zero. 
4910-5060  This  is  the  subroutine  which  selects  the  longitude  endpoints  of  the  rhumb  line 
approximation  to  a  great  circle  between  a  starting  position  with  latitude  PX 
and  longitude  LX  and  a  final  position  with  latitude  PY  and  longitude  LY.  It 
then  computes  and  prints  the  heading  and  cumulative  distance  along  each  leg. 
4910-4950  Using  the  increment  size,  IN  and  IC,  compute  the  number  of  interior  points, 
NP,  on  the  rhumb  line  approximation  exclusive  of  the  initial  and  final  points 
(the  number  of  interior  legs  is  NP  -  1.  Redefine  IC  if  it  is  longer  than  the 
difference  of  the  initial  and  final  longitudes. 
4960-4980  Find  the  heading  and  distance  on  the  initial  leg.  Print  out  first  line  of  output. 
4990  If  there  are  no  internal  legs,  go  to  5060. 

5000  Initialize  the  counter  IC.  If  there  is  only  one  leg  remaining,  go  to  5040. 
5010-5030  Loop  to  compute  print  each  internal  leg. 

5040-5050  Compute  the  headings  and  distances  for  the  final  leg.  If  M%  is  one,  return  to 
compute  the  limiting  latitude  leg. 
5060  Print  the  final  summary  line  of  output. 
51 10-5200  Subroutine  to  print  one  line  of  output.  NL  counts  the  number  of  lines  of  output. 

If  20  lines  are  output,  a  new  screen  is  started. 
5220-5260  Subroutine  to  print  screen  heading. 

5500  Return  to  BASIC. 
6000-6030    TRESS  ANY  KEY  TO  CONTINUE9  subroutine. 
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IV.  NAVEPHMx  Almanac  and  Ephemeria  Program. 

A.  Introduction.  The  basis  for  NAVBPHM  are  the  equations  of  VanFlandern  and  Pulkki- 
nen  (VFP)  [Ref.  10].  These  equations  can  be  used  to  compute  the  mean  heliocentric 
positions  of  the  sun  and  planets  and  the  mean  geocentric  position  of  the  moon  for  the 
mean  ecliptic  and  equinox  of  date.  The  authors  claim  their  formulas  to  be  of  low-precision 
(1')  and  valid  for  any  epoch  within  300  years  of  the  present.  When  corrected  for  nutation 
and  aberration  the  accuracy  of  their  formulas,  at  least  for  the  sun,  moon  and  navigational 
planets,  appears  to  be  much  better,  i.e.  0.2',  at  least. 

The  formulas  of  VFP  are  used  to  compute  heliocentric  spherical  ecliptic  coordinates 
for  any  specified  ephemeris  time,  ET.  These  coordinates  are  longitude  A,  latitude  /?  and 
radius  vector  r.  These  coordinates  must  be  converted  to  rectangular  coordinates  i,  y ,  and 
z  using  the  standard  transformation 

x  =  r  cos /9  cos  A, 

y  =  rcos/?sinA,  and  (30) 

z  —  rsin/?. 

To  obtain  the  geocentric  i,  y,  and  z  coordinates  of  the  planets,  subtract  the  x,  y,  and  z 
coordinates  of  the  sun  from  the  x,  y,  and  z  coordinates  of  the  planets,  respectively. 

The  Julian  day  number,  required  in  many  calculations,  is  obtained  using  the  equation 
on  page  B2  of  Ref.  2.  The  ephemeris  time  is  obtained  from  the  universal  time,  UT,  from 
the  equation  ET  =  UT  +  AT.  The  factor  AT  is  obtained  by  astronomical  observation 
only.  The  formula  used  here,  AT  =  81.94T  -  15,  was  obtained  by  least  squares  from  the 
1980.5  through  1985.5  values  given  on  pages  B5  and  K9  of  The  Astronomical  Almanac 
1985  [Ref.  11],  where  T  is  the  number  of  Julian  centuries  elapsed  from  1900  January  0, 
12*  ET.  Although  a  plot  of  AT  as  a  function  of  time  is  linear  for  1980.5  through  1985.5, 
this  should  be  checked  with  each  new  edition  of  The  Astronomical  Almanac.  An  accurate 
value  of  AT  affects  only  the  computation  of  the  moon's  position.  Errors  of  as  much  as 
9.6'  in  AT  will  affect  the  computation  the  moon's  position  by  at  most  0.1',  consequently 
it  will  not  be  necessary  to  change  FMPL  distribution  tapes  yearly. 

The  heliocentric  in- plane  velocity  components,  iw  and  yw,  of  the  planets,  required  for 
the  aberration  correction,  can  be  computed  from  the  formulas  [Ref.  12,  pg.  85]: 

xv  =  —aeEsmEj  and  . 

yw  =  aE\/l-e2co8Ei 
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where  E  is  the  eccentric  anomaly,  E  is  the  time  derivative  of  E,  a  is  the  semimajor  axis  and 
c  is  the  eccentricity.  Unfortunately,  to  find  E  one  mnst  solve  Kepler's  equation  iteratively, 
which  is  a  slow  process.  With  an  error  of  no  more  than  2%  for  the  navigational  planets, 
one  can  use  the  mean  anomaly  M  in  place  of  E  in  Equs.  (31).  The  in-plane,  velocity 
components  are  thus  approximated  by 


xv  =  —aeEsmM,  and 
yv  =  aEy/l-ePcosM. 

From  Equ.  (3.76)  of  Ref.  12,  it  can  be  shown  that 


(32) 
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where  /*  is  the  reduced  mass  and  k  is  the  gravitational  constant.  The  heliocentric  ecliptic 
velocity  components,  i,  y  and  i,  can  then  be  obtained  by  the  transformation 


X 

"cosft 

—  sinft 

0" 

"1 

0 

0 

*COS  (J 

—  sinw 

0 

• 

y 

ss 

shift 

cosft 

0 

0 

cosi 

-sin* 

sinw 

C08W 

0 

9w 

z 

0 

0 

1 

0 

sint 

cost 

0 

0 

1 

0 

where  Q  is  the  longitude  of  the  ascending  node,  i  is  the  orbital  inclination  and  u  is  the 
argument  of  perihelion.  The  equations  for  computing  a,  e,  Af,  u,  i  and  0  were  obtained 
from  Escobal  [Ref.  13,  pp.  8-9].  The  geocentric  i,  y  and  z  coordinates  of  the  planets  are 
obtained  by  subtracting  the  x,  y  and  z  coordinates  of  the  sun  from  the  heliocentric  i,  y 
and  z  coordinates  of  the  planets.  Since  the  aberration  of  the  moon  is  negligible,  its  velocity 
components  are  not  computed. 

Once  the  mean  geocentric  rectangular  positions  (xmi  ym  and  zm)  and  velocities  (im, 
ym  and  im)  have  been  obtained,  the  longitude  and  latitude  for  the  mean  ecliptic  and 
equinox  of  date  can  be  obtained  by  inverting  Equs.  (30),  so  that 
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(33) 
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Next  the  nutation  of  longitude  A^  and  obliquity  Ae  are  computed  [Ref.  4,  §2C].  The 
geocentric  rectangular  positions  (it,  yt  and  zt)  for  the  true  ecliptic  and  equinox  of  date 
can  be  obtained  from  the  transformation 
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The  same  transformation  is  used  for  the  velocity  components.  The  longitude  and  latitude 
for  the  true  ecliptic  and  equinox  of  date  can  be  obtained  by  substituting  the  true  positions 
into  Equs.  (33). 

The  aberration  or  light-time  correction  converts  true  positions  into  apparent  positions. 
The  x-coordinate  transformation  is 

. 

TXt 

c 

where  r  is  the  geocentric  distance  and  c  is  the  velocity  of  light.  Similar  transformations 
are  used  for  the  y-  and  z-coordinates.  The  apparent  longitude  and  latitude  for  the  true 
ecliptic  and  equinox  of  date  can  be  obtained  by  substituting  the  apparent  positions  into 
Equs.  (33). 

The  apparent  geocentric  rectangular  positions  (xaqi  y^  and  Zaq)  for  the  true  equator 
and  equinox  of  date  are  obtained  from  the  transformation 
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The  apparent  right  ascension  a  and  declination  6  for  the  true  equator  and  equinox  of  date 
can  be  computed  from 


a  =  qatn(yo*,  Xaq),  and 
6  =  sin"1  ^  =  tan-1 


»aq 
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The  Greenwich  mean  sidereal  time,  Greenwich  apparent  sidereal  time,  Greenwich 
hour  angle  and  local  hour  angle  are  computed  from  formulas  given  on  pages  B3  and  B4 
of  Ref.  2.  The  altitude  and  azimuth  are  computed  from  Equs.  (10)  and  (11),  respectively. 
The  equation  for  refraction  is  Equ.  (3),  page  B15  of  Ref.  2.  The  equations  for  planetary 
magnitude  are  given  on  pg.  315  of  Ref.  4  with  the  correction  for  Saturn's  rings  given  on 
pages  362  to  365.  Formulas  for  the  table  on  page  365  were  obtained  using  least  squares. 
Equations  for  the  semidiameter  of  the  sun  and  moon  are  found  on  page  B16  of  Ref.  2.  The 
lunar  parallax  in  altitude  formula  is  on  page  B16  of  Ref.  2  and  the  lunar  phase  formula  is 
on  page  311  of  Ref.  4.  The  lunar  age  approximation  developed  by  the  author  is  within  ±1 
day. 
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B.  Sample  Problem.  Calculate  the  position  of  the  sun,  moon  and  navigational  plan- 
ets for  24  May  1984  at  17h58m458  zulu  time  at  latitude  38°35'24"  north  and  longitude 
76°18'00"  west.  Assume  that  the  temperature  and  pressure  are  the  default  values  of  10°  C. 
and  1010  mb.,  respectively.  The  output  shown  in  this  sample  problem  was  computed  using 
Microsoft  basioa/D  (double  precision)  on  an  IBM  PO.  Several  intermediate  results,  such 
as  the  true  position,  and  apparent  positions  for  the  true  equator  and  for  the  true  equinox 
of  date  are  printed  only  as  a  debugging  aid  for  those  wishing  to  reproduce  the  results  on 
a  computer  other  than  those  for  which  NAVEPHM  is  available. 
Input  Screen: 

YEAR  (4  DIGITS)  ?  1984 

MONTH  NUMBER  ?  5 

DAY  OF  THE  MONTH  7  24 

ZULU  TIME  (HH.MMSS)       ?  17.5845 
LAT  (DD.MMSS)     (+N/-S)?  38.3524 
LON  (DDD.MMSS)    (+W/-E)?  76.18 
TEMP  (DEC  CELSIUS)         ?  10 
PRESSURE  (MILLIBARS)     ?  1010 

PRESS  C  TO  CONTINUE 

The  input  screen  is  used  for  input  only — there  are  no  options  from  which  to  choose.  The 
temperature  and  pressure  are  used  only  for  the  atmospheric  refraction  computation. 
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1st  Output  Screen: 

SUN 

TRUE  POS,  TRUE  ECL  k  EQNX  OF  DATE  (UT) : 
63.73236355163271    0 

APP  POS,  TRUE  ECL  k  EQNX  OF  DATE  (UT) : 
63.72674538386791    0 

APP  POS,  TRUE  EQU  k  EQNX  OF  DATE  (UT) : 
4 . 11446999613666     20 . 89927908414652 

GREENWICH  HOUR  ANGLE  *  DECL  (UT) : 
90°29.0'  20*54.0' 

ALTITUDE  -    68.50534784507299  68° 30. 3' 

AZIMUTH    -    218.6611413863376  218° 39. V 

SD  -  15.8' 
REFRACTION  -     .4' 

PRESS  C  TO  CONTINUE 

The  tine  and  apparent  positions  for  the  true  ecliptic  and  equinox  of  date  printed  are  the 
ecliptic  longitude  and  latitude,  respectively,  in  decimal  degrees.  The  apparent  position 
for  the  true  equator  and  equinox  of  date  is  the  right  ascension  in  decimal  hours  and  the 
declination  in  decimal  minutes.  The  Greenwich  hour  angle  and  declination  are  given  in 
degrees,  minutes,  and  tenths  of  minute  notation.  The  altitude  and  azimuth  are  given  in 
both  decimal  degrees  and  degrees,  minutes,  and  tenths  of  minute  notation.  SD  is  the  sun's 
semi-diameter.  The  refraction  is  for  the  computed  altitude  of  the  sun. 
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2nd  Output  Screen: 
MOON 


TRUE  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
365 . 6933313649804  -4 . 991896604201667 

APP  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
365 . 6933313649804  -4 . 991896504201667 

APP  POS,  TRUE  EQU  k  EQNX  OF  DATE  (UT) : 
23 . 86924633529027  -6 . 291911880374428 

GREENWICH  HOUR  ANGLE  k  DECL  (UT) : 
154°09.7'   -  6°17.5' 


ALTITUDE  -  5.45172646204884 

5°27.1' 

AZIMUTH  -  267.4666215802323 

257°28.0' 

PHASE  .31    AGE  -  24  DAYS 

SD  -  14.8'   SD  AUG  -  14. 8* 

P  IN  A  -  54' 

REFRACTION  -  10. 2? 

PRESS  C  TO  CONTINUE 

The  output  for  the  moon  is  similar  to  that  for  the  sun.  In  addition,  the  moon's  phase, 
approximate  age,  augmented  semi-diameter  (topocentric  semi-  diameter)  and  parallax  in 
altitude  P  IN  A  are  output.  The  true  position  and  the  apparent  position  for  the  true  ecliptic 
and  equinox  of  date  for  the  moon  are  identical  because  no  aberration  correction  is  made 
for  the  moon  (see  Section  A). 
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3rd  Output  Screen: 

VENUS 

TRUE  POS,  TRUE  ECL  k  EQNX  OF  DATE  (UT) : 
57.72453552418941  -.6528408164408445 

APP  P08,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
57 . 71230545404549  - . 6632006744515710 

APP  POS,  TRUE  EQU  k  EQNX  OF  DATE  (UT) : 
3 . 706665486856183  10 . 01586622158487 

GREENWICH  HOUR  ANGLE  k  DECL  (UT) : 
96°36.0'     19°01.0/ 

ALTITUDE  -  63.67703088674409  63° 40. 6' 

AZIMUTH  -  227.7083425077327  227° 42. 4' 
MAGNITUDE  -  -3.4 
REFRACTION  -  .5' 

PRESS  C  TO  CONTINUE 

The  output  for  Venus  is  similar  to  that  for  the  sun  except  that  the  apparent  magnitude  is 
output  instead  of  the  semi-diameter. 
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4th  Output  Screen: 
MARS 

TRUE  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
226 . 1549997906566  - . 7543527638263081 

APP  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
226 . 1561434601825  - . 7542106991415879 

APP  POS,  TRUE  EQU  *  EQNX  OF  DATE  (UT) : 
14 . 89744688277315  -17 . 39602620686829 

GREENWICH  HOUR  ANGLE  k  DECL  (UT) : 
288°44.3'   -  17°23.8' 

ALTITUDE  -  -54.68434644133615      -  54°41.1' 
AZIMUTH  -  62.30703743152624       62° 18. 4' 
MAGNITUDE  -  -1.6 
REFRACTION  -  -.7' 

PRESS  C  TO  CONTINUE 

The  output  for  Mara  is  similar  to  that  for  Venus.  Note  that  the  altitude  of  Mara  is  negative, 
that  is,  Mara  is  below  the  horizon.  The  leads  to  the  anomolous  computation  of  a  negative 
value  for  the  refraction.  When  a  planet  is  below  the  horizon,  the  value  computed  for 
refraction  is  meaningless. 
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5th  Output  Screen: 

JUPITER 

TRUE  POS,  TRUE  ECL  k  EQNX  OF  DATE  (UT) : 
282 . 0081844949666  . 1412184643412081 

APP  POS,  TRUE  ECL  k  EQNX  OF  DATE  (UT) : 
282 . 0100408579675  . 1412662061186586 

APP  POS,  TRUE  EQU  k  EQNX  OF  DATE  (UT) : 
18 . 86941225466823  -22 . 75882790181082 

GREENWICH  HOUR  ANGLE  k  DEGL  (UT) : 
229°09.5'   -  22°45.5' 

ALTITUDE  -  -61.97025716907849      -  61° 58, 2' 
AZIMUTH  -  296.4711031451715        296°28.3' 
MAGNITUDE  -  -2.1 
REFRACTION  -  -.5' 

PRESS  C  TO  CONTINUE 

The  output  for  Jupiter  is  similar  to  that  for  Venus.  See  the  comments  for  Mars  regarding 
a  negative  refraction  value. 
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6th  Output  Screen: 

SATURN 

TRUE  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
221 . 5236884913894  2 . 580801361264901 

APP  POS,  TRUE  ECL  *  EQNX  OF  DATE  (UT) : 
221 . 5271566155949  2 . 58091611295858 

APP  POS,  TRUE  EQU  k  EQNX  OF  DATE  (UT) : 
14 . 66067525585222  -12 . 83630963737889 

GREENWICH  HOUR  ANGLE  *  DEGL  (UT) : 
292°  17. 4'   -  12°50.2' 

ALTITUDE  -  -49.04279937791162      -  49°02.6' 
AZIMUTH  -   60.93733251413396        60° 56. 2' 
MAGNITUDE  -  .4 
REFRACTION  -  -.8' 

PRESS  C  TO  CONTINUE 

The  output  for  Saturn  is  similar  to  that  for  Venus.  See  the  comments  for  Mars  regarding 
a  negative  refraction  value. 
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C.  Program  Lieting— NAVBPHM 

10  REM  YANFLANDERN  ft  PUUXINEN  PLANET  EPHDfERIS.  07-20-85. 
REV  08-27-85  C  0800 

12  REN  "ASTRO  EPHEN  Y" 

13  REN  "NAVEPHM.BAS" 

20  DEFDBLA-N.O-Y:DIM  A(26),B$<5) 

30  P2-2*ATN(1)  :P>P2*P2:TP«PI*PI:RD-PI/180:8D-RD/3600 

40  EK-.0172020d806:MB$-n010"iDC$— 10-:DN$-CHR|(26) 

60  DEFPNA(I)-a-INT(X))*TP:DEFFNBa)-X-TP*INT(X/TP) : 
DEFFNR(X)-INT(X*MO+ .  5)/M0 

60  DEFFNM(X)=X-M0*INT(X/M0) 

70  B$(0)-"SUN" :B$ (1)-"M0QN" :B$(2)-"VENUS" :B$(3)-"MARS" 

80  B$(4)-"  JUPITER"  :B$(5)-"SATURN" 

90  GOTO  2000 

100  : 

110  REN  2-ARG  ATAN  FCTN.  RADIANS. 

120  A2-ATN(Y/a-lE-09*a-0)))-PI*a<0)*TP*a>0)*(Y<0)  :RETURN 

130  : 

140  REN  ARCSIN  *  ARCCOS  FCTN'S 

160  AS-ATN(X/(SQR(1-X*X)-1E-09*(ABS(X)-1)))  :RETURN 

160  AOATN(SqR(l-X*X)/(X-lE-09*(X«0)))-PI*(X<0)  :RETURN 

170  : 

180  REN  HELIO  3PHER-T0-RECT  ♦  NOTION  COMPUTATION 

190  R-R/100000:CS-R*C08(B) 

200  XH-C3*C0Sa):YB-CS*SIN(L):ZH-R*8IN(B) 

210  IF  N-1THEN  RETURN 

220  MU-1:IF  NoOTHEN  MU-1+1/RM(N) 

230  FA-fX*SQR(AM(N)*MU)/R 

240  X— FA*SIN(M(N))  :Y»FA*8QR(1-EC(N)*2)*C0S(M(N))  :Z-0 

250  A—  AP(N):G08UB  280:A—  IN(N)  :G0SUB  290:A—  AN(N):GOSUB  280 

260  UH-X;VH-Y:WH-Z: RETURN 

270  : 

280  CA-COS (A) :  3A-8IN (A) :  >X*CA+Y*SA :  Y-Y*CA-X*8A : X-T :  RETURN  :  REM  Z-AXIS  ROT 

290  CA=C08(A)!3A=8IN(A):T=Y*a+Z*8A:Z=Z*CA-Y*8A:Y=T;RETURN  : REM  X-AXIS  ROT 

300  : 

310  REN  RECT-TO-SPBER 

320  GOSUB  120 

330  L-A2/RD:R>X*X+Y*Y:B-ATN(Z/8QR(R2))/RD:R-8qR(R2+Z*Z) : RETURN 

340  : 

350  REN  DECIMAL  TO  ODD  MN.F 

360  V$»"  ":IF  X<OTHEN  V$-"-":X— X 

370  X=X+1/1200:Y=INT(X):V*=V$+RIGHT$("  "+STR*(Y),3)*"°" 

380  X-600* (X-Y) : Y-INT(X) : X|»8TR$( 1000+Y) 

390  Vt=V$+MID$(X$,3,2)+"."+RICTT$(X$.l):RETURN 

400  : 

410  REN  DDD.NN88  TO  DECIMAL 

420  IX-0:FOR  Z-1T0  LENCVDtCl-MIDKVl.Z.DtIF  C$-"."THEN  IX-Z 

430  NEXT  :IF  IX-OTBEN  X-VAL(V$) : RETURN 

440  X-VAL(LEFT$(V$ ,  EC) )  :  SN-1 :  IF  X<OTHEN  SN— 8N  :X— X 

460  V$-V$+"0000":Y-VAL(MID|(V$,IX+1.2))  :Z-VAL(MID|(V$,IX*3,2)) 

460  X-SN* ( (Z/60+Y) /60+X) : RETURN 
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470  : 

480  REM  FUNDAMENTAL  ARGUMENTS 

400  Aa)-.6O6434+.O366OliO129*TS:A(2)-.374897+.O36291d47O0*TS 

500  A(3)». 250001+ . 0367481962*TS  :AC6)=0:AC7)= .  779072+ . 00273790031*  T8 

510  A(8)-.903126+.OO2737778S*TS 

520  A(  12)-. 506408+4 .450468670000003E-03*TS : A( 13)  » . 140023+ . 00446036173* TS 

530  A(14)«. 292498+. 00445040017*TS:A(15)-. 987353+. 00145575328*TS 

540  A( 16) - . 063856+ . 00145561327* TS : AC 17) -  840694+ . 00145669466*13 

550  A(18)«.080608+.00023080893*TS:A(19)-.066631+.00023080893*TS 

560  A(20)=.814704+.00023080803*T8:A(21)=.133295+.00000294371*TS 

570  AC22)- . 882987+ . 0000929437 1*TS : A(23)- . 821218+ . 00009294371 *TS 

580  A(26)-.4OO580+3.26O438OOOOOOOO1E-O5*TS 

590  A(4)-AC1)-AC7):AC5)-AC1)-AC3) 

600  FOR  Z=1T0  25:A(Z)=FNA(A(Z)):NEXT  Z:RETURN 

610  : 

620  PRINT  BICO). -REM  SUNCO) 

630  I=AC3):Y=AC13):Z=AC19) 

640  L-(6010-17*TJ)*8IN(X)+72*SIN(2*X)-7*COS(X-Z)+6*8IN(A(1)-A(7)) 

660  L«L+5*8INC4*X-8*A(16)+3*Z)-6*C08(2*(X-Y))-4*8IN(X-Y) 

660  L-L+4*C08(4*X-8*A(16)+3*Z)+3*SIN(2*(X-Y))-3*8IU(Z)-3*8IN(2*CX-Z)) 

670  L=L*SD+A(7):B=0 

680  R»100014-1875*COS(X)-14*C08(2*X) 

690  AM(0)-1.00000023:EC(0)«.016709114:M(0)-A(8) 

700  AP(0)-( (4 . 62778E-04*TJ+1 . 719175) *TJ+101 . 2208333) *RD: 

IN(0)-0:AN(0)-180*RD 

710  RM(O)-1/32030O 

720  GOSUB  190 : RR-ft : XS-XH : YS-YH : ZS-ZH : US-UH : : VS-VH : W8-VH :  RETURN 

730  PRINT  B$(1):REM  MOON(l) 

740  W*A(2) :X=A(3) : Y*A(4) : Z-AC8) 

760  Irf2640*Sm(W)-4686*8IN(y-2*Y)+2370*8IN(2*Y)+769*SIN(2*W) 

760  L»L-668*8IN(Z) -412*8IN(2*X) -212*8IN(2*  CW-Y) ) -206*8INO»-2*Y+Z) 

770  >L+192*8IN(W+2*Y)+166*8IN(2*Y-Z)+148*SIN(W-Z)-125*8IN(Y) 

780  L-L-110*8IN(W+Z)-55*8IN(2*(X-Y))-46*8IN(W+2*X)+40*SIN(W-2*X) 

790  L-L-38*8IN(W-4*Y)+36*8IN(3*W)-31*SIN(2*W-4*Y)+28*SIN(W-2*Y-Z) 

800  L=L-24*8IN(2*Y+Z)+19*8INCW-Y)+18*SIN(Y+Z)+15*8INCW+2*Y-Z) 

810  >L+14*SIN(2*(W+Y))+14*SIN(4*Y)-13*SINC3*W-2*Y) 

820  ACS«W+16*AC7)-18*A(12):L»L-ll*SINCAG)+9*C0S(AC)+4*TJ*CC08CAC)+8nCAC)) 

830  L-L+10*8IN(2*W-Z)+9*8IN(W-2*X-2*Y)-9*8IN(2*W-2*Y+Z) 

840  L=L-8*8INCW+Y)+8*8INC2*Y-2*Z)-8*SINC2*W+Z)-7*8INC2*Z) 

860  L»L-7*SIN(W-2*Y+2*Z)+7*8IN(A(5))-6*SINCW-2*X+2*Y) 

860  L-L-6*8IN(2*X+2*Y) -4*8IN(W-4*Y+Z) -4*8IN(2*W+2*X) 

870  L-L+3* (SIN(W-3*Y) -8INCW+2*Y+Z) -SIN(2*W-4*Y+Z)+SINCW-2*Z) 

+8IN(W-2*Y-2*Z)) 

880  L-L+2*(SINCW+4*Y)-SINC2*W-2*Y-Z)-SIN(2*X-2*Y+Z)) 

890  L-L+2*CSIN(4*W)+8IN(4*Y-Z)+SIN(2*W-Y)) 

900  B-18461*8IN(X)+1010*8IN(W+X)+1000*8IN(W-X)-624*SIN(X-2*Y) 

910  B=B-199*SIN(W-X-2*Y)-167*SIN(W+X-2*Y)+117*8IN(X+2*Y) 

920  >B+62*8INC2*W+X)+33*SINCW-X+2*Y)+32*8IN(2*W-X)-30*8IN(X-2*Y+Z) 

930  B«B-16*SIN(2*W+X-2*Y)+16*SIN(W+X+2*Y)+12*SINCX-2*Y-Z) 

940  B-B-9*SIN(W-X-2*Y+Z) -8*SIN(X+A(5) )+8*8IN(X+2*Y-Z) 

960  B-B+7*(-8IN(W+X-2*Y+Z)+8IN(W+X-Z) -SDKW+X-4*Y) ) 

960  B«B+6*(-SIN(X+Z)-SINC3*X)+8IN(W-X-Z)) 

970  B-B+5*C-8INCX+Y)-SINCW+X+Z)-8INCW-X+Z)+8IN(X-Z)+SIN(X-Y)) 
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080  B-B-K4*(SIN(3*W*X)-8INa-4*Y))i'3*(-SIN(W-X-4*Y)+8IN(W-3*X)) 

MO  >B*2*<-SIN(2*W-X-4*Y)-8Ill(3*X-2*Y)-*8IN(2*W-X*2*Y)) 

1000  B-B*2*(SIN(W-X*2*Y-Z)*SIN(2*W-X-2*Y)+8IN(3*W-X)) 

1010  R-6O36298-32ra6*C0S(W)-B79O4*C0S(W-2*Y)-46367*C0S(2*Y) 

1020  R-R-89O4*C0S (2*W) +3866*C0S(2*W-2*Y) -3237*C03(2*Y-Z) 

1030  R-R-2688*COS(W*2*Y)-23S8*COSaf-2*Y+Z)-2030*COS(,*-Z) 

1040  HFR*171fl*C08(Y)i»1671*C08(W+Z)+1247*C08(W-2*X)+704*COS(Z) 

1050  R-R+529*C08(2*Y+Z)-624*C0S(W-4*Y)+308*CQS(W-2*Y-Z)  -366*C0S(3*W) 

1060  Rrt-296*C0S(2*W-4*Y)-263*C0S(Y+Z)*249*CTS(3*W-2*Y)-221*C0S(W+2*Y-Z) 

1070  RfR+185*C08(2*X-2*Y)-161*C08(2*Y-2*Z)+147*C08(W+2*X-2*Y)-142*C08(4*Y) 

1080  R-R+130*C08(2*W-2*Y+Z)-118*C08(W-4*Y+Z)-116*C08(2*W+2*Y) 

1000  R-R-110*C0S(2*W-Z) 

1100  >L*SD+A(1):B-B*8D:R-R/234S4.3:G0SUB  180  :RP-R:  RETURN 

1110  PRINT  R$(2):REM  VENU8(2) 

1120  X-A(8):Y-A(13):Z-A(14) 

1130  L=(2814-20*TJ)*8IN(Y)-181*8IN(2*Z)*12*SIN(2*Y)-10*C08(2*(X-Y)) 

1140  L»Lt7*C0S(3*(X-Y)) 

1150  B-1221B*SIN(Z)*83*(8IN(Y+Z)+8IN(Y-Z)) 

1160  R-72335-493*C0S(Y) 

1170  >L*SD+A(12) :  B-B*SD 

1180  AM(2)-. 7233316000000001 :EC(2)-6.773041000000001E-03:M(2)-A(13) 

1190  AP(2)*((-.001386380*TJ+. 50818611) *TJ+64. 3841881) *RD:IN(2)=3.304«*RD 

1200  AN(2)-( ( . 0004 1 *TJ+ . 80085) *TJ+75 . 77064722000001) *RD : RN(2)-408523 . 5 

1210  GOSUB  190 :RP«R: RETURN 

1220  PRINT  B$(3):REM  MAR8(3) 

1230  X-AC8) : Y-AC16) : Z-A(17) : Q-A(IO) 

1240  L-(38461*37*TJ)*SIN(Y)*(2238*4*TJ)*8IN(2*Y)+181*SIN(3*Y)-52*SIN(2*Z) 

1250  L=L-22*C08(Y-2*q)-19*8IN(Y-Q)+17*C08(Y-q)+17*8IN(4*Y)-16*C08(2*(Y-q)) 

1260  L-L-H3*C08(X-2*Y)-10*8IN(Y-2*Z)-10*SIN(Y*2*Z)*7*C08(X-Y)-7*COS(2*X-3*Y) 

1270  L«L-B.*8IN(A(13)-3*Y)-5*8IN(X-Y)-B*SIN(X-2*Y)-4*C08(2*X-4*Y)*4*C08(q) 

1280  L-L+3*C0S(A(13)-3*Y)*3*8IN(2*(Y-q)) 

1200  B-6603*SINCZ)*622*8IN(Y-Z)+615*8IN(Y+Z)+64*SIN(2*Y*Z) 

1300  R-153031 - 14 170*CflS (Y) -660*C08 (2*Y) -47+C0S (3* Y) 

1310  L=L*SD+A(15);B=B*8D 

1320  AM(3)-1. 52368830: EC(3)-. 00340488700000002 :M(3)-A(16) 

1330  AP(3)-((1.3125E-04*TJ+1. 06076667) *TJ+285. 4317610000001) *R0: 

IN(3)-1.8497*RD 

1340  AN(3)-(  (-1 .  389E-06*TJ* .  7700916700000001)  *TJ*48 .  78644167)  *RD : 

RM(3) -3098710 

1350  GOSUB  100 :RP=R: RETURN 

1360  PRINT  B$(4):REM  JUPITER(4) 

1370  X-A(18):Y-A(19):Z»A(22) 

1380  L»19934*8IN(Y)*5023*TJ+2611+1093*C08(2*Y-5*Z)+601*SIN(2*Y) 

1390  L-L-479*SIN(2*Y-6*Z)-186*SIN(2*Y-2*Z)*137*SIN(3*Y-5*Z)-131*SIN(Y-2*Z) 

1400  L-L*79*C08(Y-Z)-76*C08(2*Y-2*Z)-74*TJ*C08(Y)+68*TJ*SIN(Y) 

♦66*C0S(2*Y-3*Z) 

1410  >L+63*C08(3*Y-5*Z)+53*C08(Y-5*Z)+49*8IN(2*Y-3*Z)-43*TJ*8IN(2*Y-5*Z) 

1420  L=L-37*C0S(Y)+25*SH(2*X)+25*8IN(3*Y)-23*SIN(Y-5*Z)-19*TJ*C0S(2*Y-5*Z) 

1430  L-L+17*C08(2*Y-4*Z)+17*C0S(3*Y-3*Z)-14*SIN(Y-Z)-13*8IN(3*Y-4*Z) 

-9*C0S(2*X) 

1440  L-L*9*C0S(Z)-9*SIN(Z)-9*SIN(3*Y-2*Z)*9*SIN(4*Y-5*Z)+9 

*8IN(2*Y-6*Z+3*A(2S)) 

1460  >L-8*C08(4*Y-10*Z)f7*C08(3*Y-4*Z)-7*C08(Y-3*Z)-7*SIN(4*Y-10*Z) 
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1460  L-L-7*8IN(Y-3*Z)*6*C0S(4*Y-6*Z)-6*SIlf(3*Y-3*Z)*6*C0S(2*Z) 

-4*8IN(4*Y-4*Z) 

1470  >L-4*C0S(3*Z)*4*C0S(2*Y-Z)M*C0S(3*Y-2*Z)-4*TJ*C0S(2*Y)*3*TJ*SIN(2*Y) 

1480  L-L+3*C08(5*Z)+3*C08(5*Y-10*Z)+3*8IN(2*Z)-2*8IN(2*X-Y)+2*8I]I(2*X+Y) 

1400  L-L-2*TJ*8IN(3*Y-5*Z)-2*TJ*8IN(Y-S*Z) 

1500  B-- 4602*C08(Y)*250*SIN(Y)+227-227*COS(2*Y)*30*TJ*8INCY)*21*TJ*C08(Y) 

1610  B-B*16*8IN(3*Y-5*Z)-13*8IN(Y-5*Z)-12*C0S(3*Y)*12*8IN(2*Y) 

1620  B»B+7*C0S(3*Y-6*Z)-5*C08(Y-5*Z) 

1630  R-620883-25122*<aS(Y)-604*C08<2*Y)*260*C08<2*(Y-Z))-170*COS(3*Y-5*Z) 

1640  >R-lO6*Sm2*(Y-Z))-01*TJ*SIN(Y)-84*TJ*COS<Y)+60*SIN(2*Y-3*Z) 

1660  a-a-«7*8IN(Y-5*Z)+66*SIN(3*Y-5*Z)+63*SIN(Y-Z)-51*C0S(2*Y-3*Z) 

1660  R=R-46*SINtt)-20*COS(Y-5*Z)+27*COS(Y-2*Z)-22*COS(3*Y)-21*SIll<2*Y-5*Z) 

1570  L»L*8D+X:B=B*SD 

1580  AM(4)«6.2O2661:EC(4)-.O484O42510:M(4)-A(10) 

1500  AP(4)=< (7 . 0406E-04*TJ+ . 6804316700000001) *TJ+273 . 2775417) *RD : 

IN(4)»1.3031*RD 

1600  AN(4)»( (3. 52222E-04*TJ+1 . 01063) *TJ+00. 44338611) *RD:RM(4M047. 365 

1610  GOSUB  100 :RP-R: RETURN 

1620  PRINT  B$(5);REM  SATURN(B) 

1630  W=A(10):X=A(21):Y=A(22):Z=A(25) 

1640  L-23045*SIN(Y) ♦6O14*TJ-2680*COS (2*W-5*Y) *2S07*1 177*SIN(2*W-5*Y) 

1660  L»L-826*C08(2*W-4*Y)+802*8IN(2*Y)+42S*8IN(W-2*Y)-220*TJ*C0S(Y) 

1660  L-L-153*C08(2*W-6*Y)-142*TJ*8IH(Y)-114*C08(Y)"H01*TJ*SIN(2*W-5*Y) 

1670  L-L-70*C08(2*X)+67*8IN(2*X)+66*8IN(2*W-6*Y)+60*TJ*COS(2*W-5*Y) 

1680  L»L*41*SIN(W-3*Y)*30*8IN(3*Y)+31*SIN(W-Y)*31*8IN(2*(W-Y)) 

1600  L»L-20*COS(2*W-3*Y)-28*SIN(2*W-6*Y-*3*Z)*28*COS(W-3*Y) 

1700  L-L*22*TJ*8IN(2*W-4*Y)-22*SIN(Y-3*Z)+20*SIN(2*W-3*Y) 

1710  >L*2O*COS(4*W-1O*YM0*COS(2*Y-3*Z)+10*8IN(4*W-1O*Y) 

1720  >L-17*TJ*C08(2*Y)-16*C08(Y-3*Z)-12*8IN(2*W-4*Y)+12*C08(W) 

1730  L=L-12*8IN(2*(Y-Z))-11*TJ*SIN(2*Y)-11*C0S(2*W-7*Y) 

1740  L"L+1O*8IN(2*Y-3*Z)+1O*COS(2*(W-Y))+0*SIN(4*W-0*Y) 

1750  L-L-8*SIN(Y-2*Z)-8*a)S(2*X+Y>8*C0S(2*X-Y)*8*C08(Y-Z) 

1760  >L-8*8IN(2*X-YM*SIN(2*X+Y)-7*C0S(W-2*Y)-7*C0S(2*Y) 

1770  L=L-6*TJ*8IN(4*W-10*Y)+6*TJ*C08(4*W-10*Y)+6*TJ*SIN(2*W-6*Y) 

1780  >L-5*8IN(3*W-7*Y) -5*C0S(3* (W-Y) ) -6*C0S(2* (Y-Z) ) ♦5*SIN(3*W-4*Y) 

1700  L»L+6*SIN(2*W-7*Y)-*4*SIN(3*af-Y))+4*SIN(3*W-5*Y) 

1800  L-L+4*TJ*<»S(W-2*Y)+3*TJ*C0S(2*W-4*Y)+3*C08(2*W-6*Y+3*Z) 

1810  L-L-3*TJ*SIN(2*X)*3*TJ*C0S(2*W-6*Y)-3*TJ*C0S(2*X) 

1820  L»L+3*C08(3*W-7*Y)*3*C08(4*W-0*Y)+3*8IN(3*W-6*Y) 

1830  L=L+3*8IN(2*W-Y)+3*SIN(W-4*Y)+2*C0S(3*(Y-Z)) 

1840  >L+2*TJ*8IN(W-2*Y)+2*SIN(4*Y)-2*C0S(3*W-4*Y) 

I860  L=L-2*C0S(2*W-Y)-2*SIN(2*W-7*Y+3*Z)*2*C0S(W-4*Y)*2*C0S(4*W-11*Y) 

1860  >L-2*SIN(Y-Z) 

1870  B=8297*8IN(Y)-3348*C08(Y)+462*8IN(2*Y)-188*C08(2*Y)+186 

1880  B-B+70*TJ*CQS(Y)-71*CQS(2*W-4*Y)+46*SIN(2*W-6*Y) 

1800  B-B-45*C08(2*W-6*Y)  f  20*8111(3*  Y)-20*C08(2*W-3*Y) 

1000  B-B*18*TJ*8IN(Y)-14*C08(2*W-5*Y)-11*C08(3*Y)-10*TJ 

1010  B-B^*SIN(W-3*Y)+8*SIN<W-Y)-6*8IN(2*W-3*Y)+5*8IN(2*W-7*Y) 

1020  B-B-5*C0S(2*W-7*Y)*4*SIN(2*W-5*Y) -4*TJ*SIN(2*Y) 

1030  B=B-3*C0S(W-Y)+3*C0S(W-3*Y)+3*TJ*SIN(2*W-4*Y) 

1040  B"B+3*SIN(W-2*Y)+2*SIN(4*Y)-2*C0S(2*(W-Y)) 

1060  R-065774-53252*C08(Y)-1878*8IN(2*W-4*Y)-1482*C0S(2*Y) 

1060  R-R*817*SIN(W-Y)-530*CO8(W-2*Y)-624*TJ*8IN(Y) 
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1970  R-R*349*SIN(2*W-6*Y)*347*8IN(2*W-6*Y)*328*TJ*C0S(Y) 

1980  R*R-225*SIN(Y)+149*C08(2*W-6*Y)-126*C0S(2*(W-Y)) 

1900  RfR+104*C0S(W-Y)*101*C0S(2*W-5*Y)*98*C0S(W-3*Y) 

2000  R»R-73*C0S(2*W-3*Y)-62*C08(3*Y)+42*SIN(2*Y-3*Z) 

2010  RfR+41*8IN(2*(W-Y))-40*8IN(W-3*Y)+40*C08(2*W-4*Y) 

2020  RpR-28*TJ-23*8IN(W)+20*8IN(2*W-7*Y) 

2030  L»L*SD*X:B-B*8D:L9-L:B9-B 

2040  AM(5)-9. 554747000000001: EC (B)-. 05554609270000001 :M(B)-A(22) 

2060  AP(5)=( CO. 73542E-04*TJ+1.086220«0)*TJ+338. 3077722) *RD:IN(5)=2.4386*RD 

2060  AN(5)-((-l. 52181E-04*TJ+. 8731961400000001) *TJ+112. 7903889) *RD: 

RM(5)-3498.5 

2070  GOSUB  190 :BP«R: RETURN 

2080  : 

2090  CLS  : INPUT  aYEAR(4  DIGITS)  ■;! 

2110  INPUT  "MONTH  NUMBER  a;M 

2130  INPUT  "DAY  OF  THE  MONTH  M 

2150  INPUT  "ZULU  TIME  (HH.MMSS)  a;V|:GOSUB  420:UT=X 

2170  INPUT  "LAT  (DD.MMS8)  (+N/-S)  a;V$: GOSUB  420:L>X 

2190  INPUT  "LON  (D0D.MM8S)  (+W/-E)  a;V$:GOSUB  420:LG=X 

2210  INPUT  "TEMP  (DEC  CELSIUS)  »;DC 

2230  INPUT  "PRES8URE  (MILLIBARS)  ";MB 

2240  GOSUB  2720 

2260  JD-367*K-INT(7*a+INT((M*9)/12))/4)*INT(275*M/9)*I+1721013.5 

2260  CLS 

2270  TJ«(JD-2415020)/36626:DT-81 .94*TJ-16 : TS-JD-2451546+DT/3600/24 

2280  T0=T8/36625 : TS=TS*UT/24 : TJ=TS/36626+l 

2290  TM$-« (UT) : ■ : IF  DT-OTHEN  TMt»« (TDT) : ■ 

2300  GOSUB  490 

2310  GM-(2 . 58622E-05*T0+2400 . 051336) *T0+6 . 697374560000001 +1 .  002737909*UT 

2320  M0=24:GM=FNM(GM):GS=GM-.0OO2O*SIN(A(5)) 

2330  DL-<-17.23-.02*TJ)*SIN(A(5))-1.27*SIN(2*A(7))*.21*SIN(2*A(5)) 

2340  DL=DL-.2*8IN(2*A(1)):DL=DL/3600*RD:REM  DL=NUTATION  OF  LONGITUDE 

2360  OE»((.00181*TJ-. 0069) *TJ-46. 846) *TJ+84428. 26+9. 21*C0S(A(5)) 

2360  0E=(0E+.552*CO8(2*A(7)))*8D 

2370  FOR  N=OT0  5: ON  N+1G0SUB  620,730,1110,1220,1360,1620 

2380  IF  N=O0R  N=1THEN  X=XH:Y=YH:Z=ZH:IF  N=1THEN  U=0:V=0:W=0:GOTO  2410 

2390  IF  N-OTHEN  U^JH:Y-VH:W-WH:XS-X:YS-Y:ZS-Z:US«U:VS-V:WS-W:GOTO  2410 

2400  X-XS+XH : Y-YS+YH : Z-ZS+ZH : U-US+UH : V-VS+VH : W-WS+WH 

2410  A=-DL:G08UB  280 

2420  PRINT  DN$;aTRUE  POS,  TRUE  ECL  ft  EQNX  OF  DATE  a;TM$: GOSUB  320:PRINT  L,B 

2430  GD=R:L1=L:B1=B:IF  N=OTHEN  LO=L 

2440  0173 .  142 :  ROR/C :  REM  VELOCITY  OF  LIGHT 

2450  X=X-RC*U : Y=Y-RC*V : Z=Z-RC*W 

2460  PRINT  :PRINT  aAPP  POS,  TRUE  ECL  ft  EQNX  OF  DATE  a;TM$: 

GOSUB  320: PRINT  L;B 

2470  A=-OE: GOSUB  290: GOSUB  320 : LP*L : BP=B : IF  N=OTHEN  LS=L;BS=B 

2480  L-L/15 

2490  PRINT  :PRINT  aAPP  P08,  TRUE  EQU  ft  EQNX  OF  DATE  a;TM|:PRINT  L;B 

2500  M0-360:GH-FNM(15*(GS-D) 

2510  PRINT  : PRINT  "GREENWICH  HOUR  ANGLE  ft  DECL  a;TM$ 

2520  X»GH:GOSUB  360 :L$-V$:X-B: GOSUB  360:PRINT  L$+"'";SPC(3);  V$+-/a 

2530  LH-(GH-LG)*RD:CL-COS(LH) :SL-8IN(LH) :BD-B*RD:SB-SIN(BD) :CB-COS(BD) 
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2540  LD-LT*RD:a-C0S(LD);SI-8IN(LD):I-ai*SB-K3*CB*CL:G0SUB  150: 

A9-AS:AL»AS/RD 

2550  Y--CB*SL:X-SB*CX-CB*8X*CL: GOSUB  120:AZ»A2/RD 

2560  PRINT  :PRINT  "ALTITUDE  »  ";AL; :X=AL:G0SUB  360:PRINT  SPC(5);V$+"'" 

2565  PRINT  "AZIMUTH  -  ";AZ; :X-AZ:G0SUB  360: PRINT  SFC(5);V$+*'« 

2570  ON  NGOSUB  2780.2780,2780,2860,2870 

2580  IF  N»1THEN  PRINT  "PHASE" ;X;"  AGE  »";AG;"  DAYS 

2500  IF  N>1THEN  M0=10. PRINT  "MAGNITUDE  =  ";FNR(MG) 

2600  IF  N^OTHEN  S«16.994/R:M0»10:PRINT  "SD  »";STR*(FNR(S))+"'" 

2610  IF  NOITHEN  2660 

2620  0*23464. 8*R: 8*936. 75/D:M0»10: PRINT  "SD  *";STR$(FNR(S)  )+"'"; 

2630  8-8* (1+SIN(A9)/D) .PRINT  ■  SD  AUG  ""jSTRKFNRO))*"'" 

2640  X-C0S(A9)/D: GOSUB  150:PA-AS/RD*6O 

2650  PRINT  "P  IN  A  ■  ";STR|(FNR(PA))*"'" 

2660  RF-1/TAN((AL+7.31/(AL+44.4))*RD) 

2670  RF-RF* ( (MB- 80 ) /930) / ( 1 + . 00008* (RF+39) * (DC- 10) ) 

2680  RF-RF-.06*SIN((14.7*RF+13)*RD):MO-10:RF-FNR(RF) 

2600  PRINT  "REFRACTION  »  ";STR$(RF)+a/" 

2700  GOSUB  2720:CLS  :NEXT  N:G0T0  2000 

2710  : 

2720  PRINT  :  PRINT  "PRES8  C  TO  CONTINUE" 

2730  d-INXEY$  :IF  C$-""THEN  2730 

2740  IF  C$-"C"THEN  RETURN 

2760  GOTO  2730 
2770 

2780  CD-8IN(B3*RD) *SIN (BP*RD) +C0S(B8*RD) *COS(BP*RD) *COS( (LS-LP) *RD) 

2790  X*CD: GOSUB  160 :Y*RR*SIN (AC) :X*GD-RR*CD: GOSUB  120: PA-ABB (A2/RD) 

2800  ON  NGOTQ  2810,2840,2850 

2810  K-(l-KBS(A2))/2:SN-SIN((L0-Ll)*RD):X-8QRa):G0SUB  150 

2820  AO-20.53*AS/PI:IF  3N>0THEN  AG-29.53-AG 

2830  M0-100:K*FNR(K)  :M0»1  :AG*FNR(AG)  :  RETURN 

2840  MG-(PA*PA*4 .  247E-07+ . 01322) *PA+2 . 1715*L0G(GD*RP) -4 : RETURN 

2860  MG*  01486*PA+2.1715*L0G(GD*RP)-1.3:RETURN 

2860  HG=2. 1715 *LOG(GD*RP)- 8. 93: RETURN 

2870  BQ»(168. 1176+1. 394091*TJ)*RD:II-(28.0743-.0127991*TJ)*RD 

2880  NN«( ( . 242202*TJ+3 . 98599) *TJ+126 . 3629) *  RD 

2890  J>(  ( .  0171866*  T  J- .  4549142)  *TJ+6 .  912936)  *RD 

2900  SQ-(- . 239992*TJ-2 . 731279) *TJ+42 . 92039 

2910  CB-C0S(B9) :  Y-SIN(II)*SIN(B9)-fC0S(II)*CB*SIN(L9-B0) 

2920  X-CB*C08(L9-B0) : GOSUB  120:UP-A2/RD 

2930  D-BP*RD:CB*C08(D) :SB*8IN(D) :8J-8IN(JJ) :CJ*C08(JJ) 

2940  AG=LP*RD-NN:SN=SIN(AG) 

2960  Y«8 J*8B+CJ*CB*8N : X*CB*C08 (AG) .GOSUB  120 : UU-A2/RD 

2970  BB*SJ*CB*SN-CJ*SB 

2980  MG*2 .  1715*L0G(GD*RP)  -8 .  68+ .  044*AB8(UP+S0-UU)  -2 .  6*ABS(BB)  +1 .  25*BB*BB 

2990  RETURN 


46 


D.  Program  Annotation— NAVEPHEM 

Line(s)  Usage 

30  P2  =  ir/2,  PI  =  t,  TP  =  2x,  RD  =  degree-to-radian  conversion. 
40  EK  as  Sun '8  gravitational  constant.  MB$  =  atmospheric  pressure  in  millibar's 
(1010  is  the  default  value),  DCS  =  temperature  in  degrees  Celsius  (10  is  the 
default  value). 
50  FNA  converts  angles  in  revolutions  to  degrees  between  0  and  2t,  FNR  rounds 

to  the  nearest  M0. 
60  FNM  is  the  X  modulus  M0  function. 
70-80  B$  is  an  array  of  planet  names. 
90  GOTO  2090  to  commence  execution. 
120  Two  argument  arctangent  function  with  output  interval  of  0  to  2x.  The  -1E-9 

must  be  replaced  by  +1E-9  in  the  PC-1500A. 
150  The  arcsine  function  can  be  replaced  inline  by  the  ASN  function  in  the  PC- 

1500A. 
160  The  arccosine  function  can  be  replaced  inline  by  the  ACS  function  in  the  PC- 
1500A. 
190*200  TJnscale  the  radius  vector  by  a  factor  of  1E5.  Convert  spherical  coordinates  to 
rectangular  coordinates. 
210  Return  if  the  body  is  the  moon. 
220-260  MU  is  the  reduced  mass  of  the  body.  Compute  the  approximate  heliocentric 
velocity  vector  of  the  body. 
280  Z-axis  rotation. 
290  X-axis  rotation. 
320-330  Convert  rectangular  coordinates  to  spherical  coordinates.  L  =  longitude,  5  = 

latitude  and  R  =  distance. 
360-390  Convert  decimal  degrees  to  degrees,  minutes  and  tenth  minute  notation. 
420-460  Convert  DDD.MMSS  or  HH.MMSS  format  to  decimal. 
490-600  Compute  the  fundamental  planetary  arguments  [Ref.  10] . 
620-720  Compute  the  geocentric  spherical  position  and  velocity  of  the  Sun.  Convert  to 

geocentric  rectangular. 
720-1100  Compute  the  geocentric  spherical  and  rectangular  coordinates  of  the  Moon. 
1110-1210  Compute  the  heliocentric  spherical  position  and  velocity  of  Venus.  Convert  to 

heliocentric  rectangular. 
1220-1350  Compute  the  heliocentric  spherical  position  and  velocity  of  Mars.  Convert  to 

heliocentric  rectangular. 
1360-1610  Compute  the  heliocentric  spherical  position  and  velocity  of  Jupiter.  Convert 

to  heliocentric  rectangular. 
1620-2070  Compute  the  heliocentric  spherical  position  and  velocity  of  Saturn.  Convert  to 
heliocentric  rectangular.  Saturn's  longitude  and  latitude  are  saved  as  L9  and 
B9  on  line  2030. 
2090-2240  Input  prompts  should  be  designed  so  that  default  values  or  previously  entered 
values  are  displayed.  Pressing  the  return  key  will  input  a  displayed  value.  K 
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any  value  is  changed,  then  the  entire  value  must  be  keyed  in.  This  feature 

had  to  be  removed  on  the  TRS-SO  Model-4  version  since  it  was  difficult  to 

implement. 
2250  Computation  of  the  Julian  Day  Number. 
2270  TJ  =  number  of  Julian  centuries  from  noon,  January  0,  1900  to  midnight  of 

input  date.  DT  =  AT  correction  factor.  [NOTE:  This  equation  should  be 

examined  for  accuracy  yearly — see  text).  TS  =  number  of  Julian  days  from 

noon,  January  1,  2000  to  0000  hours  Universal  Time  of  the  input  date. 
2280  TO  =  number  of  Julian  centuries  from  noon,  January  1,  2000  to  0000  hours 

Universal  Time  (UT)  of  the  input  date.  TS  =  number  of  Julian  days  from 

noon,  January  1,  2000  to  the  input  date  and  time.  TJ  =  number  of  Julian 

centuries  from  noon,  January  0,  1900  to  the  input  date  and  time. 
2290  If  the  AT  correction  is  set  to  zero,  then  the  positions  are  referenced  to  TDT 

(terrestial  dynamic  time)  rather  than  to  UT. 
2300  Compute  fundamental  arguments. 
2310  GM  =  Greenwich  Mean  Sidereal  Time. 
2320  GS  =  Greenwich  Apparent  Sidereal  Time. 
2330-2340  DL  =  nutation  of  longitude. 
2350-2360  OE  =  obliquity  of  the  ecliptic  corrected  for  nutation. 

2370  FOR  loop  to  cycle  through  Sun,  Moon  and  planets.  N  =  body  number.  The 

loop  ends  on  line  2700. 
2380  X,  Y  and  Z  are  the  geocentric  coordinates  of  the  Sun  or  the  Moon.  Set  the 

velocity  components  of  the  Moon  equal  to  zero. 
2390  Save  the  position  (XS,  YS  k  ZS)  and  velocity  (US,  VS  &  WS)  components  of 

the  Sun. 
2400  Compute  the  geocentric  position  and  velocity  components  of  the  Nth  planet. 
2410  Correct  position  for  nutation. 

2420  Display  true  position  for  the  true  ecliptic  and  equinox  of  date. 
2430  Save  true  distance  GD,  longitude  LI,  latitude  Bl  and  solar  longitude  L0. 
2440-2450  Correct  position  for  aberration  (light  time). 

2460  Display  apparent  position  for  the  true  ecliptic  and  equinox  of  date. 

2470  Convert  ecliptic  coordinates  to  equatorial  coordinate.  Save  the  right  ascension 

LP  and  declination  BP  of  the  body.   For  the  Sun,  save  right  ascension  and 

declination  as  LS  and  BS. 
2480  Convert  right  ascension  from  degrees  to  hours. 
2490  Display  apparent  position  for  the  true  equator  and  equinox  of  date. 
2500  GH  =  Greenwich  Hour  angle. 
2510-2520  Display  Greenwich  Hour  angle  and  declination. 
2530-2565  Compute  and  display  the  altitude  and  azimuth.  The  altitude  is  saved  as  A9  on 

line  2540. 
2570  Subroutine  call  for  N  =  1  through  5. 
2580  Display  the  Moon's  phase  and  age. 
2590  Display  the  planet's  magnitude. 
2600  Compute  and  display  the  Sun's  semidiameter  [Ref.  2,  pg.  B16]. 
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2610 

2620 

2630 

2640-2650 

2660-2690 

2700 
2720-2760 
2780-2790 

2800 
2810 
2820 
2830 
2840 
2850 
2860 
2870-2980 

2870-2900 


2910-2920 

2930-2970 

2980 


Transfer  if  the  body  is  not  the  Moon. 

Compute  and  display  the  Moon's  semidiameter  [Ref.  2,  pg.  B16]. 

Compute  and  display  the  Moon's  augmented  semidiameter  [Ref.  2,  pg.  B16). 

Compute  and  display  the  Moon's  parallax  in  altitude  [Ref.  2,  pg.  B16). 

Compute  and  display  the  refraction  correction  at  the  body's  altitude  [Ref.  2, 

pp.  B14-B15). 

Continue  prompt.  End  of  FOR-NEXT  loop  on  body  number  N. 
Continue  prompt  query. 
Compute  the  phase  angle  PA.  CD  is  the  cosine  of  the  elongation  [Ref.  4, 

pg.  312). 

Transfer  for  the  Moon,  Venus  or  Mars. 
Compute  the  Moon's  phase  angle  K  [Ref.  4,  pg.  311]. 
AG  is  an  approximation  to  the  Moon's  age  developed  by  the  author. 
Round  the  Moon's  age  to  the  nearest  day  and  return. 
MG  is  the  magnitude  of  Venus  [Ref.  4,  pg.  314] . 
MG  is  the  magnitude  of  Mars  [Ref.  4,  pg.  314]. 
MG  is  the  magnitude  of  Jupiter  [Ref.  4,  pg.  314]. 
Compute  the  magnitude  MG  of  Saturn.  This  computation  is  complicated  by 

the  varying  aspect  of  Saturn's  rings.  Details  follow. 
Computation  of  the  ring's  ascending  node  BO,  inclination  II,  right  ascension 

NN,  inclination  to  the  mean  equator  JJ  and  arc  SO  from  NN  to  BO-  These 

equations  result  from  curve  fitting  the  table  in  Reference  4,  page  365. 
Computation  of  UP  =  U'  [Ref.  4,  pg.  364]. 
Computation  of  UU  =  U  and  BB  =  sin  B  [Ref.  4,  pg.  345]. 
MG  is  the  magnitude  of  Saturn  [Ref.  4,  pg.  314]. 
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APPENDIX:  The  QATN  Function 

This  routine  is  the  standard  arctangent  function  corrected  for  quadrant.  The  quadrant  arc- 
tangent function  is  occasionally  implemented  as  the  ATAN2  function,  the  ANGLE  function 
or  the  Rectangular-to-Polar  function. 

Entering  variables  are  the  x-  and  y-coordinates,  X  and  Y.  The  exiting  variable  is 
the  angle  6,  where  —  t  <  8  <  r.  Use  of  the  quadrant  arctangent  function  is  denoted  by 
e  =  qatn(7,JC). 

1.  If  X  ^  0,  go  to  step  4. 

2.  Sete  =  (*/2)*sgn(y). 

3.  Go  to  step  8. 

4.  Set  0  =  arctan(y/X). 

5.  If  X  >  0,  go  to  step  8. 

6.  Sete  =  e  +  **sgn(y). 

7.  if  y  =  o,  set  e  =  *■. 

8.  Return. 

Note: 

If  y  >  0  then  sgn(y)  =  +1. 
If  y  =  0  then  sgn(y)  =  0. 
If  y  <  0  then  sgn(y)  =  -1. 

Users  of  Microsoft  BASIC  can  simplfy  the  qatn  function  significantly  by  using  the  code 
given  below.  To  return  an  angle  of  9  (designated  by  A  in  the  code)  in  the  range  of  (-ir,  jt), 
use: 

PI  -  4*ATN(1):     TP  =■  PI  +  PI:  EPS  -  1E-33 

A  -  ATN(Y/(X-EPS*(X-0)))  -  PI*(X<0)*(SGN(Y)   -  (Y=0)) 

To  return  a  value  of  A  in  the  range  of  (0, 2ir),  use: 

PI  -  4*ATN(i):     TP  -  PI  +  PI:  EPS  -  1E-33 

A  -  ATN(Y/(X-EPS*(X»0)))  -  PI*(X<0)  +  TP*(X  >=  0)*(Y<0) 
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